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1.  Overview:  Brain  Structure-Function  Couplings 


1.1  Objective 

This  overview  aims  to  capture  the  general  direction  of  this  evolving  program  as  developed  by  the 
collaborators  on  the  initiative. 

1.2  Program  Overview 

Every  individual  has  a  unique  set  of  talents  and  skills,  and  these  individual  differences  manifest 
in  the  speed  and  accuracy  of  behavioral  responses  across  tasks,  the  strategies  employed  to  solve 
problems,  and  the  reactions  to  challenging  experiences,  to  name  just  a  few.  Individuals  have 
unique  capabilities,  and  researchers  have  attempted  to  explain  these  differences  in  human 
behavior  for  centuries.  Many  disciplines  have  relied  on  behavior  alone,  but  the  field  of 
neuroscience  capitalizes  on  continual  advancements  in  brain  imaging  methodologies, 
computational  approaches,  and  modeling  techniques  to  better  understand  and  predict  the 
interaction  between  the  brain  and  behavior.  By  studying  neural  processing  mechanisms  in  the 
brain,  neuroscientists  hope  to  unravel  how  the  brain  enables  the  mind. 

Recent  advancements  in  neuroimaging  technologies  and  computational  analysis  capabilities  have 
enabled  new  approaches  for  understanding  the  basic  principles  of  brain  organization  and 
function.  The  research  initiative  on  brain  structure-function  couplings  capitalizes  on  some  of 
these  advancements  to  examine  how  differences  in  brain  structure  and  brain  function  can  be 
linked  to  differences  in  human  behavior  (figure  1). 


Figure  1 .  The  initiative  investigates  how  couplings  between  individual  variation  in  brain  structure  and  brain 
function  can  be  linked  to  individual  differences  in  human  behavior. 

The  brain  consists  of  both  grey  matter,  primarily  neuronal  cell  bodies  and  glial  cells,  and  white 
matter,  myelinated  axons  that  transmit  electrochemical  signals  between  grey  matter  regions.  The 
grey  matter  is  heterogeneous.  Different  regions  of  the  brain,  often  referred  to  as  Brodmann’s 
areas  after  the  scientist  who  pioneered  and  popularized  a  particular  parcellation  scheme,  have 
different  cellular  properties  and  structure,  and  the  field  of  cytoarchitectonics  examines  how  grey 
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matter  is  structured.  A  complementary  field  of  myeloarchitectonics  examines  how  axons  are 
structured.  Although  grey  matter  regions,  such  as  the  aforementioned  Brodmann  areas,  are 
largely  preserved  between  people,  these  regions  will  differ  in  size,  shape,  location,  and 
cytoarchitecture  across  individuals.  Likewise,  major  fiber  tract  bundles  of  axons  that  connect 
particular  grey  matter  regions  are  preserved  across  people,  but  the  tracts  within  the  bundle  vary 
in  the  total  number,  projection  path,  termination  location,  and  myeloarchitecture  of  the  tracts 
between  individuals.  These  numerous  sources  of  individual  structural  variation  hold  potential  for 
explaining  individual  differences  in  neuronal  activity  within  a  grey  matter  region  and  neural 
communication  between  regions.  Critically,  recent  advancements  in  structural  neuroimaging 
technologies,  known  as  diffusion  weighted  imaging  techniques,  provide  a  new  capability  to 
image  white  matter  fiber  tracts  noninvasively  in  human  volunteers.  These  new  white  matter 
imaging  techniques  can  be  combined  with  existing  robust,  structural  imaging  measures  for  grey 
matter.  Within  the  past  decade,  the  neuroscience  field,  including  this  research  initiative,  is  just 
beginning  to  examine  approaches  that  best  utilize  these  structural  imaging  techniques  in 
conjunction  with  functional  imaging  data  of  brain  activity  and  measurements  of  human  behavior. 

Similar  to  structural  variation,  there  are  also  numerous  sources  of  functional  variation  between 
people.  The  electrochemical  activity  within  a  region  has  spatial  topologies  and  temporal 
dependencies  dependent  on  the  cytoarchitecture  of  the  grey  matter  itself.  These  variations  reflect 
concentrations  of  neurotransmitters,  organization  of  synapses,  refractory  periods  for  these 
various  electrochemical  processes,  etc.  Likewise,  the  transmission  time  and  projection  path  of  a 
neuronal  signal  depends  on  the  myeloarchitecture  of  the  fiber  tracts,  as  well  as  lower  level 
electrochemical  properties  of  the  tracts  themselves.  These  structural  constraints  of  functional 
activity  are  expected  across  these  spatial  scales  within  a  physical  system,  but  the  current  imaging 
technologies  cannot  easily  measure  the  available  resources  and  dynamic  structure  at  each  of 
these  levels.  Consequently,  our  research  efforts  emphasize  functional  parameters  and  structure 
that  can  be  imaged  noninvasively  in  human  subjects,  and  in  particular,  this  research  initiative 
also  exploits  recent  developments  in  functional  connectivity  analysis  techniques.  Traditionally, 
functional  analysis  methods  have  focused  on  individual  differences  in  localization  of  brain 
activity.  For  example,  studies  may  identify  whether  a  particular  grey  matter  region  is  activated 
across  people,  or  whether  the  timing  of  a  neural  response  within  a  region  is  the  same  between 
different  tasks.  Functional  connectivity  approaches  emphasize  the  global  dynamics  and  topology 
of  brain  activity.  While  some  functional  connectivity  research  has  looked  at  coupled  activity 
between  pairs  of  brain  regions,  recent  innovative  approaches  have  started  examining  the 
properties  of  larger  networks  and  computing  connectivity  across  multiple  regions 
simultaneously.  These  newly  developed  methods  emphasize  a  systems-level  or  network-level 
description  of  communication  across  the  brain,  providing  an  opportunity  to  leverage  insights 
from  network  science. 

Combined,  the  new  imaging  methods  for  structural  connections  and  the  novel  methodologies  for 
understanding  global  functional  connectivity  provide  a  promising  avenue  to  examine  the  link 
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between  structural  variation  and  functional  variation  across  individuals.  These  developments 
hold  promise  to  identify  individual  variations  that  can  be  linked  to  differences  in  human 
behavior.  Importantly,  this  initiative  emphasizes  the  development  of  imaging  measurements  and 
analysis  methods  that  capture  statistically  meaningful  differences  between  individuals.  If  the 
predictive  capability  between  structure,  function,  and  behavior  is  realized,  brain  metrics  can  be 
used  in  a  host  of  applications  ranging  from  neurotechnologies  that  optimize  Soldier-system 
performance,  to  training  technologies  that  decrease  time  to  train,  to  assessments  of  whether 
injured  Soldiers  should  return  to  battle.  In  short,  this  initiative  addresses  a  fundamental  question 
in  neuroscience  about  the  predictive  coupling  between  variations  in  brain  structure  and  brain 
function,  and  if  it  is  successful,  it  will  have  broad-based  application  to  the  Army. 


1.2.1  Near-term  Goals 


For  the  near-term  goals  of  this  initiative  in  Brain  Structure-Function  Couplings,  we  feel  it  is 
important  to  focus  on  particular  scales  of  nervous  system  organization.  Like  other  multiscale 
modeling  domains,  the  nervous  system  has  important  features  at  multiple  spatial  scales  as  well  as 
multiple  temporal  scales  within  a  given  spatial  scale.  As  visualized  in  figure  2,  the  spatial  scales 
of  brain  structure  range  from  genes  and  small  proteins  up  through  single  neurons  to  local  circuits 
within  small  areas  of  the  cortex,  all  the  way  up  to  distributed  networks  of  brain  regions  that 
enable  global,  whole-brain  interactions. 


Spatial  Scales  of  Organization 


Brain  Structure 


network  function 

functional 

connectivity  global  activation 

patterns 


synchronized  activation 


spike  trains 


bursts 


action  potential 


neurotransmitter  binding 
(chemical) 


voltage  change 
(electrical) 


depolarization  hyperpolarization 

membrane  resistance 


transcription 


translation 


Figure  2.  The  nervous  system  can  be  described  at  multiple  spatial  scales  for  both  brain  structure  (left)  and  brain 
function  (right). 
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Likewise,  brain  function  can  also  be  measured  and  described  at  these  multiple  spatial  scales. 
Recent  modeling  efforts  like  the  Blue  Brain  project  are  using  supercomputers  to  replicate  the 
structure  and  function  of  local  circuits  known  as  cortical  columns  from  thousands  of  experiments 
in  the  rat  cortex  (EFPL,  2011).  The  majority  of  the  data  is  collected  from  single  neurons  within  a 
cortical  column  that  averages  a  total  of  10,000  neurons,  and  computational  algorithms  are  then 
built  for  the  column  and  tested  to  see  how  well  the  model  replicates  the  empirical  data.  Using 
this  approach,  the  modeling  effort  intends  to  gradually  build  a  replication  of  the  entire  rat  cortex. 
Our  Director’s  Strategic  Initiative  (DSI)  effort  is  not  aimed  at  building  a  replication  of  whole- 
brain  activity  in  humans;  instead,  our  research  efforts  focus  on  the  circuit  and  distributed  system 
scales,  examining  if  network-level  brain  structural  organization  and/or  network-level  brain 
functional  activation  patterns  are  predictive  of  task-relevant  brain  states.  These  scales  were 
largely  selected  based  on  current  neuroimaging  modalities  that  measure  global  brain  activity  and 
ongoing  innovations  for  imaging  structural  networks  and  computing  functional  connectivity.  To 
this  end,  we  examine  the  circuit  and  distributed  system  scale  of  nervous  system  organization  in 
two  complementary  efforts:  an  empirical  effort  and  a  modeling  effort.  Critically,  these  two 
research  efforts,  and  the  two  studied  spatial  scales,  must  intersect.  Our  group  has  chosen  to 
leverage  the  Connectome  approach  (Bullmore  and  Spoms,  2009)  where  graph  theory  serves  as 
the  framework  for  both  structural  and  functional  network  descriptions. 

Using  a  graph  theoretic  framework,  a  common  scheme  is  used  to  generate  a  structural  brain 
graph  and  a  functional  brain  graph  to  capture  an  individual’s  network  properties.  Brain  regions 
are  considered  nodes  on  the  graph,  and  connections  between  those  regions  are  the  edges.  For 
either  type  of  graph,  the  nodes  of  the  graph  must  first  be  determined  by  dividing  the  grey  matter 
of  the  brain  into  smaller  regions,  a  process  that  relies  on  imaging  an  individual’s  brain  anatomy 
using  a  magnetic  resonance  imaging  (MRI)  scanner.  Brain  regions  used  as  graph  nodes  could  be 
very  coarse  and  segment  the  cortex  into  the  four  major  cortical  lobes  (occipital,  temporal, 
parietal,  and  frontal),  smaller  grey  matter  regions  could  be  defined  using  a  standard  atlas  (e.g., 
Harvard-Oxford  atlas),  task-specific  regions  could  be  identified  by  a  functional  dataset,  or 
regions  could  be  defined  by  numerous  other  parcellation  methods.  Once  graph  nodes  are  defined, 
the  edges  can  be  separately  defined  for  the  structural  graph  and  the  functional  graph.  The  edges 
for  an  individual’s  structural  graph  are  defined  by  the  white  matter  fiber  tract  connections  among 
the  brain  regions  used  as  graph  nodes.  An  individual’s  fiber  tracts  are  collected  using  a  diffusion- 
weighted  imaging  (DWI)  scan  on  a  MRI  scanner,  and  the  location  and  direction  of  the  white 
matter  tracts  are  inferred  from  the  constrained  movement  of  water  molecules  (see  Hagmann  et 
al.,  2006  for  a  review).  The  edges  for  a  functional  graph  are  defined  by  a  statistical  coupling 
between  the  recorded  brain  activity  across  multiple  brain  regions.  Consequently,  functional 
graphs  can  be  computed  from  time  series  of  an  individual’s  brain  activity  collected  using  any  one 
of  several  different  neuroimaging  techniques  such  as  electroencephalography  (EEG),  functional 
magnetic  resonance  imaging  (fMRI),  magnetoencephalography  (MEG),  etc.  The  end  result  is  a 
graph  with  connections  between  nodes,  providing  a  common  framework  for  both  structure  and 
function  within  individuals.  These  graphs  serve  as  the  basis  and  the  intersection  for  all  of  our 
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empirical  and  modeling  efforts,  as  discussed  in  more  detail  below.  In  addition,  these  graphs 
provide  a  mechanism  to  incorporate  temporal  scales  within  our  spatial  scale  of  nervous  system 
organization,  as  connection  strengths  between  nodes  can  be  altered  to  reflect 
leaming/experience,  brain  damage,  etc.  Thus,  the  adopted  graph  theoretic  framework  lends  itself 
to  multiscale  modeling  efforts  across  both  spatial  and  temporal  scales. 

As  described  in  1.2.1,  our  empirical  research  efforts  aim  to  identify  both  structural  and  functional 
variations  that  are  statistically  meaningful  indices  of  individual  function  and  behavior.  The  first 
is  a  crude,  but  common,  approach  that  directly  examines  the  predictive  relationship  between 
differences  in  a  structural  parameter,  such  as  fractional  anisotropy  (FA),  and  differences  in  task 
performance  (Bengtsson  et  ah,  2005;  Cascio  et  ah,  2007;  Madden  et  al,  2004;  Scholz  et  al., 

2009).  The  work  also  examines  relationships  between  functional  analysis  measures,  both 
traditional  and  new  connectivity  measures,  and  differences  in  behavior.  Thus,  our  work  extends 
what  structural  and  functional  parameters  are  explored  for  predictive  relationships  with 
behavioral  performance.  Our  second  approach  explores  predictive  links  between  structural 
parameters  and  differences  in  functional  activation  patterns,  including  the  amplitude  of  the  brain 
activity,  latency  of  the  activity,  topography  of  the  activation  pattern,  etc.  Again,  our  research 
examines  novel  pairings  of  structural  parameters  and  functional  activity  measures  as  well  as  their 
interactions  with  behavioral  performance.  Our  third  approach  capitalizes  on  the  adopted  graph 
theory  framework  for  the  project.  Numerous  graph  measures  can  be  computed  that  capture 
features  of  a  given  network,  such  as  its  modularity,  centrality,  small  worldliness,  etc.  (see  Bassett 
et  al.,  2011  for  a  review  of  graph  measures).  Our  group’s  working  hypothesis  is  that  this 
network-level  description  holds  promise  for  identifying  predictive  couplings  between  individual 
variations  in  structure  and  function  and  their  interactions  with  behavioral  performance.  Suppose 
the  modularity  of  structural  networks  correlate  with  the  ability  of  particular  functional 
connectivity  measures  to  capture  a  task-relevant  brain  activity  pattern.  Neurotechnologies  can 
then  be  designed  that  match  the  most  sensitive  analysis  method  for  the  functional  data  based  on 
the  structural  properties  of  a  Soldier’s  brain.  This  third  approach  looks  for  predictive 
relationships  within  the  properties  of  the  structural  and  functional  graphs  themselves.  Across  all 
three  approaches,  the  empirical  efforts  identify  promising  network  parameters  and/or 
connectivity  patterns  that  capture  meaningful  individual  variation  and  uses  these  parameters  to 
help  guide  the  modeling  efforts. 

This  research  initiative  contains  two  modeling  efforts:  one  is  tied  closely  with  the  empirical 
efforts  described  previously  and  focuses  on  modeling  electrochemical  processing  in  the  brain, 
while  the  second  emphasizes  the  biomechanical  response  properties  of  brain  tissue  and  hopes  to 
bridge  tissue  damage  to  electrochemical  processing  disruptions/changes.  Both  efforts  embody 
graph  theory  within  their  model.  The  electrochemical  modeling  effort  simulates  interconnected 
functional  nodes  that  produce  electrical  activity  at  frequencies  and  temporal  patterns  controlled 
through  parameterization  of  the  model.  This  effort  builds  from  David  and  Friston  (2003)  work  to 
simulate  the  functional  activity  of  small  patches  of  cortex  within  each  node.  Many  different 
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simulated  networks  can  be  implemented  that  have  different  numbers  of  nodes,  different  activity 
profiles  within  the  nodes,  different  connectivity  patterns  among  the  nodes,  communication 
between  nodes  at  various  timeframes,  etc.  Thus,  this  electrochemical  modeling  effort  provides  a 
mechanism  to  explore  how  different  structural  connections  influence  regional  activity,  functional 
connectivity,  and  global  network  properties.  It  also  provides  a  testbed  to  examine  measures  of 
network  properties.  This  modeling  effort  and  the  empirical  research  efforts  are  intended  to  refine 
and  augment  one  another,  as  one  is  used  to  constrain  or  inform  the  other. 

The  other  modeling  effort,  the  biomechanical  effort,  focuses  only  on  brain  structure.  Here,  the 
graph  nodes  represent  brain  tissue  with  information  about  the  white  matter  fiber  tracts 
incorporated  into  the  finite  element  model  to  better  model  brain  tissue  for  simulation  work. 
Research  here  is  focused  on  developing  a  damage  law  that  prescribes  the  brain  tissue  loading  and 
response  properties,  based  on  empirical  studies,  for  particular  stress  and  strain  rates  at  various 
orientations  and  durations.  Simulations  can  then  be  run  that  predict  damage  to  brain  tissue,  and 
this  damage  will  then  be  transferred  to  a  network  framework  so  the  model  can  produce  a 
simulated  structural  graph  that  is  altered  to  reflect  the  simulated  tissue  damage.  For  example,  the 
connection  strength  between  nodes  may  be  decreased  when  the  damage  law  hits  a  certain 
threshold  or  a  node  of  the  structural  graph  may  be  deleted  at  a  different  threshold  value.  By 
adopting  the  graph  framework  for  representing  simulated  damage  to  the  structural  network,  the 
biomechanical  effort  can  link  back  to  the  electrochemical  modeling  and  empirical  research.  The 
“damaged  network”  is  another  graph  with  different  connections  among  nodes,  and  the 
electrochemical  modeling  effort  provides  a  mechanism  to  explore  how  the  functional  activity 
changes  due  to  these  changes  in  the  underlying  network  connections.  This  discussion  of  the 
modeling  efforts  highlights  the  importance  of  considering  temporal  scales  within  the  study  of 
structure-function  couplings.  Certainly  in  traumatic  brain  events,  the  network  properties  are 
changing  as  the  injury  dynamics  unfold,  but  temporal  scale  is  equally  critical  in  healthy 
individuals  and  any  attempt  to  link  structural  graphs  and  functional  graphs  has  to  acknowledge 
the  substantial  differences  between  their  inherent  timescales. 

The  near-term  effort  focuses  on  structure-function  couplings  based  on  a  short-term  snapshot  of 
structural  and  functional  graphs  to  establish  the  foundational  principles;  however,  the  long-term 
goal  of  the  program  is  to  develop  time-evolving  predictive  models  of  structural  variations  that 
influence  time-evolving  models  of  global  brain  function  and  dynamic  changes  in  cognitive 
performance.  Both  structural  and  functional  connections  change  on  varying  timescales. 

Structural  adaptations  are  inherently  slow.  Structural  changes  in  professional  concert  pianists 
correlate  with  the  number  of  adolescent  practicing  hours  (Bengtsson  et  al.,  2005),  but  structural 
changes  have  also  been  found  after  a  few  days  of  intensive  training  to  juggle  (Scholz  et  al., 

2009).  An  individual’s  white  matter  pathways  also  reflect  body  mass  index  (Xu  et  al.,  2011). 
However,  structural  changes  can  occur  on  very  short  timescales  in  conjunction  with  blast- 
induced  or  blunt  force  traumatic  brain  injuries  (Taber  et  al.,  2008).  Likewise,  functional 
adaptations  can  also  reflect  changes  over  different  timescales.  The  brain  constantly  shifts 
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between  stable  brain  states  while  performing  a  task.  These  changes  can  occur  as  quickly  as 
200  ms.  Global  patterns  of  functional  activity  have  been  correlated  with  predicted  learning 
success  on  a  new  motor  task  across  days  (Bassett  et  al,  2011),  but  they  also  reflect  long-term 
states,  such  as  particular  neural  disorders  (Stam  et  al,  2007).  These  examples  highlight  the  range 
of  meaningful  timescales  already  identified  in  the  literature,  and  the  inherent  challenge  of  linking 
the  timescales  of  structural  and  functional  changes  where  the  timescales  for  measurable  changes 
are  quite  different.  Our  near-term  efforts  are  looking  at  short-term  scales,  establishing  scaleable 
frameworks  for  future,  time-evolving  research  efforts. 

1.2.2  Summary 

The  Brain  Structure -Function  Couplings  research  effort  addresses  a  basic  Army  challenge:  how 
can  we  predict  individual  Soldier  performance?  This  basic  research  (6.1)  effort  seeks  to  unravel 
the  coupling  between  structural  and  functional  variations  that  reflect  individual  differences  from 
genetics,  development,  experience,  training,  etc.,  and  link  these  differences  to  measurable 
differences  in  behavior. 

The  Army  benefits,  however,  do  not  depend  on  the  long-term  achievement  of  the  entire  project. 
Each  of  the  four  research  areas — electrochemical  modeling,  biomechanical  structural  changes, 
electrochemical  data  collection  and  analysis,  and  time-evolving  functional  connectivity — has 
been  scoped  to  develop  products  to  benefit  the  U.S.  Army  Research  Laboratory  (ARL)  mission 
programs  in  the  near-term  (figure  3).  The  empirical  work  capitalizes  on  recent  advancements  in 
both  structural  imaging  methods  and  global,  functional  connectivity  measures  in  order  to 
examine  novel  structural  and  functional  parameters.  These  efforts  are  collecting  new  datasets 
with  structural  and  functional  data  from  the  same  individuals  in  collaboration  with  academic 
partners  to  explore  couplings  within  Army  relevant  domains.  Identifying  novel  metrics  for 
individual  differences  will  aid  the  development  of  individual-specific  neurotechnologies.  In 
concert,  the  electrochemical  modeling  effort  examines  basic  principles  of  how  structural 
connections  effect  functional  relationships,  allowing  the  empirical  and  modeling  results  to  refine 
each  other’s  future  efforts.  In  particular,  it  currently  serves  as  a  critical  test  bed  to  evaluate  if 
functional  measures  developed  on  empirical  datasets  capture  the  expected  connectivity  patterns 
in  a  simulation  where  the  communication  dynamics  are  known.  Finally,  the  biomechanical 
modeling  effort  examines  empirical-based  approaches  for  developing  damage  laws  the 
incorporate  “injury  physics”  for  brain  tissue  and  ways  to  threshold  and  transfer  this  damage  to  a 
structural  connectome.  This  effort  increases  the  fidelity  of  simulations  and  can  enhance  Soldier 
protection  technologies.  Across  all  of  the  research  efforts,  this  DSI  in  Brain  Structure-Function 
couplings  builds  an  in-house  multidisciplinary  expertise  in  computational  brain  modeling  that  is 
essential  from  basic  needs  for  reviewing  proposals  to  sophisticated  efforts  to  identify  modeling 
advancements  that  can  be  leveraged  for  Army  problems 
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Figure  3.  An  overview  of  the  four  main  research  areas  in  the  DSI  on  Brain  Structure-Function  Couplings, 
a  brief  description  of  their  scientific  products,  and  how  the  areas  are  envisioned  to  intersect  with 
one  another. 

1.3  Report  Overview 

This  report  highlights  a  subset  of  Year  1  accomplishments  in  the  program.  Within  the  empirical 
efforts  on  electrochemical  processes  of  the  brain,  we  have  chosen  to  highlight  our  work  on 
functional  measures  this  year.  The  first  two  projects  describe  complementary  analyses  of  a 
previously  collected  EEG  dataset  where  Soldiers  used  a  demilitarized  rifle  in  a  target  shooting 
task  (Kerick  et  al.,  2007).  Both  projects  aim  to  evaluate  if  a  functional  connectivity  measure, 
Weighted  Phase  Lag  Index  (WPLI)  (Vinck  et  al.,  2011),  can  successfully  capture  brain  network 
dynamics  in  the  presence  of  real-world  movement  artifacts.  Project  1  is  an  exploratory  analysis 
of  one  study  participant,  completed  by  a  201 1  summer  research  intern,  Sandhya  Rawal.  This  first 
project  compared  the  WPLI  measure  to  a  traditional  time-frequency  analysis,  and  the  results 
indicated  that  the  new  connectivity  measure  captured  at  least  three  different  brain  activity 
patterns  not  identified  by  the  traditional  analysis  approach.  These  preliminary  results  suggest  that 
the  method  may  have  advantages  over  traditional  approaches  when  the  dataset  has  large-scale 
movement  artifacts  (ARL-TM-201 1,  2011).  In  a  complementary  effort,  Project  2  examines  two 
properties  of  WPLI:  its  resilience  to  movement  artifacts  (in  isolation  and  concurrent  with  task¬ 
relevant  brain  activity)  and  its  sensitivity  to  brain  activity  (again,  in  isolation  and  concurrent  with 
movement  artifact).  Across  the  12  subjects  analyzed,  the  measure  was  significant  for  brain- 
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related  activity,  but  it  was  not  significant  for  the  artifacts  (i.e.,  it  was  resilient  against  the  noise 
artifacts).  These  results  indicate  that  time-evolving  WPLI  holds  promise  for  understanding  brain 
function  in  complex  operational  environments  with  real-world  eye  and  body  movements, 
providing  a  critical  tool  for  the  development  of  neurotechnologies  that  could  monitor  neural  data 
in  real  time  to  augment  Soldier-system  performance. 

The  final  two  projects  described  in  this  report  showcase  the  accomplishments  under  the 
program’s  two  modeling  efforts — Project  3  describes  progress  under  the  electrochemical 
modeling  research  and  Project  4  describes  progress  under  the  biomechanical  modeling  research 
(Kraft  and  Dagro,  in  press).  In  Project  3,  a  neural  mass  model  that  generates  EEG-like  simulated 
functional  data  was  implemented.  Several  simulations  of  two-node  networks  verified  the 
oscillatory  behavior  of  the  simulated  brain  regions  (network  nodes).  The  simulated  EEG  data 
was  then  analyzed  using  functional  connectivity  measures,  including  the  time-evolving  WPLI 
measure  that  was  tested  on  the  empirical  datasets  in  Projects  1  and  2.  Since  the  underlying 
functional  activity  is  known,  the  simulation  provides  a  powerful  testbed  for  understanding  and 
interpreting  the  results  from  the  connectivity  measures.  The  first  round  of  simulation  results 
indicate  that  the  tested  measures  do  capture  the  temporal  dynamics  of  the  communication 
between  nodes,  but  even  in  two-node  networks,  the  results  indicate  that  network  structure 
strongly  influences  the  measure  results  and  suggests  relationships  between  structure  and  function 
need  to  be  explored  further.  Finally,  in  Project  4,  a  method  is  described  for  incorporating 
anisotropic  properties  of  biological  tissue  into  a  finite  element  model  computational  framework. 
The  algorithm  links  the  directionality  information  from  empirically  collected  structural  data, 
namely  the  white  matter  fiber  tracts  that  connect  brain  regions,  to  the  transverse  isotropic  model. 
This  work  lays  the  critical  foundation  needed  to  increase  the  fidelity  of  the  brain  tissue  response 
in  simulations  of  tissue  damage  due  to  blast  or  blunt  forces. 

Collectively,  these  four  projects  capture  the  research  enhancements  derived  through  cross- 
Directorate  collaboration.  The  simulated  EEG  data  can  be  used  to  validate  and  better  understand 
functional  measures  that  are  tested  on  empirical  datasets.  The  imaged  brain  structure  is 
incorporated  into  a  finite  element  model  to  increase  the  fidelity  of  the  simulated  tissue  response 
properties.  Critically,  in  our  first  year,  each  of  the  four  efforts  has  successfully  incorporated  the 
Connectome-based,  graph  theoretic  framework  to  ensure  advancements  between  areas  can  be 
integrated  as  the  research  develops  in  the  years  to  come. 

1.4  Year  1  Transitions  (Oct  2010-Sept  2011) 

The  following  have  been  generated: 

•  Refereed  Journal  Papers: 

1.  Lau,  T.  M.;  Gwin,  J.  T.;  McDowell,  K.;  Ferris,  D.  P.  Weighted  Phase  Lag  Index 
Stability  as  an  Artifact  Resistant  Measure  to  Detect  Cognitive  EEG  Activity  During 
Locomotion.  Journal  of  NeuroEngineering  and  Rehabilitation,  submitted. 
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Conference  Presentations: 


1.  Vettel,  J;  Bassett,  D;  Kraft,  R;  Grafton,  S.  Physics-based  Models  of  Brain  Structure 
Connectivity  Informed  by  Diffusion-weighted  Imaging.  Army  Science  Conference, 

Nov  2010. 

2.  Vettel,  J;  Bassett,  D;  Grafton,  S.  Reproducibility  of  Diffusion-weighted  Imaging 
Across  Multiple  Scanning  Sessions.  Society  for  Neuroscience,  Nov  2010. 

•  Government  Reports: 

1 .  Kraft,  R;  Dagro,  A.  Design  and  Implementation  of  a  Numerical  Technique  to  Inform 
Anisotropic  Hyperelastic  Finite  Element  Models  using  Diffusion  Weighted  Imaging ; 
(Project  4  in  this  report)  ARL-TR-5796;  U.S.  Army  Research  Laboratory:  Aberdeen 
Proving  Ground,  MD,  in  press. 

2.  Rawal,  S.  Weighted  Phase  Lag  Index  (WPLI)  as  a  Method  for  Identifying  Task-Related 
Functional  Networks  in  Electroencephalography  (EEG)  Recordings  during  a  Shooting 
Task ;  (Project  1  in  this  report);  ARL-TM-201 1;  U.S.  Army  Research  Laboratory: 
Aberdeen  Proving  Ground,  MD. 

3.  Amrein,  B.  E.;  Cassenti,  D.  N.;  Cosenzo,  K.  A.;  Krausman,  A.  S.;  LaFiandra,  M.  E.; 
Lance,  B.  J.;  Lockett,  J.  F.;  McDowell,  K.  G.;  Oie,  K.  S.;  Rice-Berg,  V.  J.;  Samms,  C. 
L.;  Vettel,  J.M.  (co-editors,  alphabetical  order).  Human  Research  and  Engineering 
Directorate:  Major  Laboratory  Programs:  Current  Thrust  Areas  and  Recent  Research', 
ARL-SR-213;  U.S.  Army  Research  Laboratory:  Aberdeen  Proving  Ground,  MD,  2010. 

•  Collaborations/Extemal  Funding: 

1 .  Jean  Vettel  (Human  Reearch  and  Engineering  Directorate  [HRED])  with  Walt 
Schneider  and  Don  Krieger  at  University  of  Pittsburgh  (Collaborative  Technology 
Alliances  [CTA]  seedling  funds  secured  for  data  collection  and  joint  analysis). 

2.  Vettel  (HRED)  with  Mike  Miller,  Barry  Giesbrecht,  and  Scott  Grafton  at  Institute  for 
Collaborative  Biotechnologies  (ICB)/University  of  California  at  Santa  Barbara  (UCSB) 
(CTA  seedlings  funds  secured  for  data  collection  and  joint  analysis  as  well  as  funds  set 
aside  for  a  joint  postdoc  to  be  hired  in  FY12). 

3.  Vettel/Kaleb  McDowell  (HRED)  funded  Paul  Sajda  at  Columbia  University  (task  order 
funds  secured  for  data  collection  and  ARL  analysis). 

4.  Reuben  Kraft  (Weapons  and  Materials  Research  Directorate  [WMRD])  sponsored 
United  States  Military  Academy  (USMA)  Davies  Fellow  Project,  “Influence  of  Brain 
Structure  in  NeuoTrauma  Protection”. 

5.  Kraft  (WMRD)  with  Barclay  Morrison  at  Columbia  University. 
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6.  Kraft  (WMRD)  with  Rika  Wright  at  Johns  Hopkins  University. 

7.  Kraft  (WMRD)  with  University  of  Pennsylvania  neural  network  Multidisciplinary 
University  Research  Initiative  (MURI). 

8.  Kraft  (WMRD)  with  Institute  for  Soldier  Nanotechnologies. 

9.  Manuel  Vindiola,  Peter  Chung,  and  Brian  Henz  (CISD)  and  Vettel  (HRED)  with 
Suvranu  De  and  Alex  Herzog  at  Rensselaer  Polytechnic  Institute. 

•  Personnel: 

1 .  Vindiola,  a  post-doc  researcher  from  Johns  Hopkins  University  hired  as  CISD  postdoc 
to  lead  electrochemical  modeling  effort. 

2.  Several  student  summer  interns  (Herzog  -  HRED/CISD,  Rawal  -  HRED,  Samantha 
Luo  -  HRED,  Philip  McKee  -  WMRD,  etc.). 

3.  Mentoring  Science  and  Math  Academy  student  on  project  “Modeling  Damage  in  the  C. 
elegans  Structural  Connectome”  (WMRD,  Kraft). 

•  Invited  Cross-Service  Briefs: 

1.  McDowell,  Medical  Research  and  Materiel  Command  (MRMC),  Feb  2011:  “ARL’s 
Neuroscience  Strategic  Research  Initiative:  Non-Injury.” 

2.  Chung,  White  House,  Office  of  Science  and  Technology  Policy  (OSTP),  March  2011: 
“Materials  for  Physical  and  Trauma  Protection  Deep  Dive.” 

3.  McDowell  and  Vettel,  Senate  Arms  Services  Committee,  Professional  Staff  Member, 
April  2011:  “ARL  Neuroscience”  and  “Brain  Structure-Function  Couplings.” 

4.  Kraft,  ARO  MURI,  July  2011:  “Computational  Brain  Biomechanics  and  Connections 
with  the  Human  Connectome.” 

5.  McDowell  and  Vettel,  Army  Laboratory  Assessment  Group,  July  2011:  “ARL 
Neuroscience”  and  “Brain  Structure-Function  Couplings.” 

•  Awards: 

1.  Kraft:  2010  Presidential  Early  Career  Award  for  Scientists  and  Engineers. 

2.  Rawal  (Mentor:  Vettel):  2nd  place  undergraduate  at  ARL  Summer  Internship 
Symposium  (c.f.  ARL-TM-201 1). 
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2.  Project  1:  WPLI  as  a  Method  for  Identifying  Task-Related  Functional 
Networks  in  EEG  Recordings  during  a  Shooting  Task 


2.1  Authors 

The  author  on  this  report  is  Sandhya  Rawal  (ARL-TM-201 1,  2001).  Significant  contributors  on 
the  project  were  Jean  Vettel  (mentor),  Samantha  Luo,  and  Kaleb  McDowell  from  ARL/HRED. 

2.2  Introduction 

EEG  records  electrical  activity  from  the  brain  by  placing  electrodes  along  the  scalp,  and  the 
study  of  this  activity  through  EEG  is  a  prevalent  methodology  used  in  neuroscience.  EEG  has 
several  advantages  over  other  prevalent  methodologies,  in  that  it  can  be  used  outside  of  the 
laboratory  and  has  a  high  temporal  resolution.  EEG  data  have  been  analyzed  at  individual 
electrode  sites  to  identify  patterns  of  brain  activity  when  an  event  occurs.  However,  the  brain 
relies  on  communication  between  brain  regions,  so  developing  analysis  methods  to  investigate 
brain  activity  between  regions  may  provide  additional  information  over  analyzing  electrical 
activity  at  single  locations.  Specifically,  identifying  the  communication  between  brain  regions 
that  occurs  during  tasks  may  provide  information  regarding  the  cognitive  processes  involved  in 
the  tasks.  The  brain  regions  and  their  connectivity  involved  in  this  communication  during  a  task 
are  referred  to  as  the  task-related  functional  network. 

One  approach  to  identify  task-related  functional  networks  is  to  investigate  the  synchronization  of 
neural  activity  between  different  brain  regions.  A  common  approach  calculates  statistical 
interdependencies  between  electrical  activity  recorded  from  pairs  of  electrodes.  Each  electrode 
can  capture  electrical  data  from  different  local  regions  of  the  brain,  so  the  interdependencies 
between  two  electrodes  may  suggest  neural  communication.  Although  neural  communication  is 
fast,  occurring  in  milliseconds,  EEG  has  a  high  temporal  resolution,  making  it  a  suitable  method 
for  identifying  communication  among  brain  regions.  However,  EEG  often  has  a  low  signal-to- 
noise  ratio  based  on  movement  and  muscle  artifacts.  Another  phenomenon  that  contributes  to  the 
low  signal-to-noise  ratio  is  often  referred  to  as  the  volume  conduction  problem,  when  two  nearby 
electrodes  tend  to  pick  up  activity  from  the  same  source,  giving  rise  to  false  correlations  between 
the  time  series. 

To  address  these  problems,  Stam  and  colleagues  (2007)  introduced  a  measure  called  Phase  Lag 
Index  (PLI)  to  address  volume  conduction  issues.  PLI  is  based  upon  the  idea  that  nonzero  phase 
lag  between  two  time  series  is  not  caused  by  volume  conduction  from  a  single  source.  Therefore, 
it  may  capture  the  communication  between  brain  regions,  as  these  communications  have  a  delay, 
or  a  nonzero  phase  lag,  due  to  axonal  transmission  properties.  The  delay  exists  because  brain 
regions  must  communicate  through  a  physical  biological  medium,  neurons,  and  thus  cannot 
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communicate  instantly.  PLI  is  less  affected  by  common  source  problems  than  the  traditional 
measures  of  functional  connectivity. 

However,  according  to  Vinck  and  colleagues  (201 1),  PLI  is  limited  by  the  discontinuity  of  phase 
lag  estimates  between  zero  and  infinitesimally  small  non-zero  phase  lags.  They  claimed  that  both 
the  sensitivity  to  noise,  including  volume-conduction,  and  the  capacity  to  detect  changes  in  phase 
synchronization  is  hindered  by  the  discontinuity  of  PLI.  To  improve  upon  PLI,  Weighted  Phase 
Lag  Index  was  introduced  by  Vinck  and  colleagues.  In  WPLI,  the  contribution  of  the  observed 
phase  leads  and  lags  is  weighted  by  the  magnitude  of  the  imaginary  component  of  the  cross¬ 
spectrum.  Vinck  and  colleagues  (2011)  demonstrated  two  advantages  of  WPLI  over  PLI,  reduced 
sensitivity  to  noise  sources  and  increased  power  to  detect  changes  in  phase  synchronization. 

The  advantage  offered  by  WPLI  of  eliminating  noise  in  EEG  data  is  particularly  useful  for  Army 
purposes  because  it  allows  for  future  studies  and  product  development  based  on  noisy  EEG  data 
that  would  be  encountered  in  real-life  situations  outside  of  the  laboratory.  Artifacts  introduced  in 
real-world  settings  seriously  exacerbate  the  low  signal-to-noise  ratio,  but  the  ability  to  identify 
cognitive  activity  within  the  noisy  data  may  allow  for  discovery  of  cognitive  processes  in  real- 
world  settings.  Identifying  nonzero  phase  lag  statistical  interdependencies  via  WPLI  may  allow 
for  the  identification  of  this  cognitive  activity. 

WPLI,  therefore,  presents  itself  as  a  methodology  to  provide  functional  connectivity  data  in 
studies  performed  in  real-world  conditions.  Accordingly,  this  project  applies  WPLI  analysis  to 
previously  collected  EEG  data  from  a  Marine  performing  a  shooting  task.  WPLI  was  computed 
on  epochs  of  data  surrounding  the  trigger  pull  event,  which  consists  of  both  brain  activity  and 
movement  artifacts  from  the  weapon  recoil.  This  project  evaluates  the  use  of  WPLI  as  a  method 
to  provide  functional  connectivity  data  that  identifies  task-related  cognitive  activity  not  identified 
using  traditional  time-frequency  analyses. 

2.3  Methods  or  Approach 

2.3.1  Participants 

In  the  original  study  (Kerick  et  al.,  2007),  18  participants  volunteered,  and  17  of  the  participants 
were  U.S.  Marines  and  1  was  a  U.S.  Army  Ranger.  All  had  completed  basic  training  and 
shooting  qualifications.  The  data  collected  from  1  participant  of  the  18  was  selected  for  analysis. 

2.3.2  EEG  Acquisition 

The  data  were  previously  collected  by  Kerick  et  al.  (2007).  The  data  were  acquired  at  500  Hz 
using  a  Neuroscan  system  with  34  active  electrodes,  including  2  horizontal  and  2  vertical  eye 
channels  and  2  mastoid  reference  channels.  These  analyses  use  the  28  non-reference  electrodes. 
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2.3.3  Task  and  Event  Extraction 

The  EEG  data  analyzed  in  this  report  were  acquired  from  the  participants  during  a  simulated 
shooting  task.  The  original  study  by  Kerick  and  colleagues  (2007)  involved  eight  shooting 
conditions.  The  conditions  included  single  and  dual  task  conditions  in  which  the  Soldier  was 
required  to  perform  only  a  shooting  task  or  perform  a  shooting  task  while  performing  an 
arithmetic  task;  with  either  short  (2-4  s)  and  long  target  (4-6  s)  exposure  and  with  all  enemy 
targets  or  mixed  enemy  and  friendly  targets  (Kerick  et  al.,  2007).  However,  this  report 
investigates  only  the  condition  that  had  a  single  task  load,  all  enemy  targets,  and  long  target 
exposure. 

The  shooting  scenarios  were  a  simulation  of  an  outdoor  shooting  range  at  Aberdeen  Proving 
Ground,  MD,  and  the  rifle  provided  to  the  participant  was  demilitarized  to  emit  a  laser  beam  for 
fire.  The  rifle  also  simulated  the  recoil  that  follows  an  actual  shot  in  order  to  imitate  the  type  of 
artifact  obtained  in  a  real-world  shooting  scenario.  A  digital  terminal  board  sent  event  markers, 
including  onset  and  offset  times  of  targets,  trigger  responses,  and  target  hits  and  misses,  from  the 
simulation  system  to  the  EEG  recording  system.  Since  all  targets  were  enemy  targets,  no 
decision  to  shoot  was  necessary.  The  participant  had  to  search  for  the  enemy  target,  orient  the 
weapon  towards  the  target,  aim,  and  pull  the  trigger. 

Thirty  trials  of  the  event  type  in  which  participants  aimed  and  shot  successfully  were  extracted 
by  epoching  the  data  around  the  event  from  -3  to  2  s,  where  0  s  corresponds  to  the  event  marker 
from  the  digital  terminal  board,  indicating  that  a  shot  was  fired  at  the  enemy  target. 

2.3.4  EEG  Pre-processing 

All  pre-processing,  including  bandpass  filtering  and  some  data  rejection,  was  performed  and 
described  by  Kerick  and  colleagues  (2007).  EEGLAB  scripts  were  used  to  automatically  remove 
bad  channels.  Independent  components  (ICs)  were  obtained  from  independent  component 
analysis  (ICA)  decomposition.  Those  with  eye  movements  or  muscle  artifacts  were  removed. 

The  channel  data  without  the  eye  movement  and  muscle  artifact  ICs  were  then  analyzed. 

2.3.5  Analysis 

WPLI  and  time-frequency  data  in  6  frequency  bands,  delta  through  gamma,  were  calculated.  The 
derived  WPLI  and  time-frequency  was  extracted  in  the  epochs  surrounding  the  trigger  pull  (-3  s 
to  2  s).  For  each  electrode  (pair)  and  frequency  band,  significance  was  determined  by  comparing 
each  time  point  (-3  to  0  s)  across  all  30  shots  to  a  baseline  determined  by  averaging  across  the 
whole  data  epoch.  The  significance  level  was  obtained  by  dividing  0.05  by  501 
(p  <  9.98003992  x  10  5)  to  account  for  the  chance  occurrence  of  a  significant  test  (Type  I  error). 
To  further  correct  for  the  Type  I  error,  at  least  five  consecutive  time  points  had  to  be  significant 
in  order  to  be  considered  a  significant  point  grouping  for  WPLI  and  at  least  25  consecutive  time 
points  were  required  for  significance  in  the  time-frequency  analysis.  The  requirement  for 
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consecutive  time  points  for  significance  is  based  on  the  estimated  time  course  of 
communications  between  brain  regions. 

2.4  Results  and  Discussion 

WPLI  and  time-frequency  measures  were  calculated  from  EEG  data  acquired  from  a  Soldier 
during  performance  of  a  shooting  task.  This  project  attempts  to  identify  task-related  functional 
networks  to  further  understand  cognitive  processes  that  occur  before  the  trigger  pull;  however, 
the  movement  artifact  caused  by  the  simulated  recoil  makes  the  signal-to-noise  ratio  low.  A 
technique  that  is  impervious  to  artifacts  of  this  type  may  allow  for  identification  of  cognitive 
processes  occurring  concurrently  with  the  movement  artifact.  In  this  project,  we  investigate 
WPLI  to  determine  if  task-related  functional  networks  can  be  identified  after  removing  the  zero 
phase  lag  statistical  interdependencies  between  pairs  of  channels.  A  traditional  time-frequency 
analysis  at  individual  channel  locations  is  used  as  a  point  of  comparison. 

The  channels  that  were  identified  as  significant  in  time-frequency  data  are  indeed  different  from 
the  ones  identified  in  WPLI.  Twenty-seven  groups  of  points  were  significant  prior  to  the  trigger 
pull  in  the  28  electrodes  and  6  frequency  bands  studied  using  time-frequency  analysis.  Fifty- 
eight  groups  of  points  were  significant  in  the  378  electrode  pairings  in  the  6  frequency  bands 
studied  using  WPLI  analysis.  In  figure  4,  the  results  are  plotted  on  a  schematic  view  of  the 
electrodes  on  the  head,  where  the  front  of  the  brain  is  located  at  the  top  of  the  graph  (Fz).  In  (a), 
all  electrodes  are  present  and  labeled  on  the  figure  in  red.  The  Fz  electrode  is  marked  by  a  larger 
red  circle.  The  significant  channels  that  were  identified  in  time-frequency  data  are  highlighted  in 
black.  Figure  4b  illustrates  the  significant  channel  pairs  identified  by  WPLI  with  black  lines 
between  the  independent  electrodes.  The  location  of  each  electrode  in  (b)  is  the  same  as  in  (a). 


Figure  4.  (a)  The  map  of  electrodes  depicted  highlights  the  electrodes  with  significance  in  black,  (b)  The  map  of 
electrodes  depicted  highlights  the  nonzero  phase  lag  statistical  interdependencies  identified  by  WPLI 
with  black  lines  between  the  interdependent  electrodes.  The  location  of  each  electrode  is  the  same  as  in  a. 
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The  differences  between  figure  4a  and  b  highlight  the  additional  information  provided  by  the 
WPLI  measure.  The  significance  of  individual  channel  locations  indicated  by  the  electrodes 
highlighted  in  black  in  figure  4a  appears  to  be  distributed  symmetrically.  However,  the 
distribution  in  the  WPLI  analysis  in  figure  4b  is  more  heavily  localized  to  the  left  hemisphere. 
Interestingly,  the  WPLI  analysis  reveals  a  left  lateralization  of  the  paired  network 
communication  among  brain  regions,  providing  an  insight  not  captured  in  the  channel  analysis. 
The  involvement  of  the  frontal  motor  cortex  and  the  parietal  cortex  suggests  that  the  cognitive 
activity  identified  may  be  motor-spatial  planning,  as  the  frontal  cortex  is  involved  in  planning 
and  in  motion  and  the  parietal  cortex  integrates  sensory  information  from  different  modalities, 
particularly  for  spatial  sense.  The  increased  interdependencies  in  the  left  hemisphere  is  consistent 
with  previous  literature,  as  left  temporal  alpha  power  is  consistently  observed  in  marksmen  few 
seconds  before  the  trigger  pull  (Kerick  et  ah). 

In  addition  to  the  lateralization,  the  WPLI  analysis  also  revealed  activity  in  the  occipital 
electrodes  not  found  in  the  channel-based  analysis.  All  occipital  electrodes  show  interaction  with 
the  F7  electrode,  suggesting  communication  between  the  visual  cortex  and  the  frontal  cortex. 

This  communication  may  reflect  that  the  shooting  task  requires  planning  in  response  to  a  visual 
target  before  the  trigger  pull.  However,  the  possibility  that  communication  between  the  occipital 
electrodes  and  the  F7  electrode  reflect  the  eye  movements  occurring  near  the  F7  electrode  must 
be  addressed  in  further  research,  as  this  report  does  not  address  the  potential  issue  of  brain 
activity  correlated  to  muscle  activity. 

Figure  4  highlights  the  promise  of  WPLI  to  reveal  network  communication  among  brain  regions. 
The  plots  show  all  pairs  of  electrodes  that  are  statistically  significant  point  groupings,  but  we 
were  also  interested  in  what  pairs  reveal  an  extended  window  of  communication  during  the  3  s 
before  the  trigger  pull.  As  a  preliminary  exploration,  we  plotted  the  four  longest  groups  of 
significant  periods,  which  are  illustrated  in  figure  5.  The  blue  line  represents  the  WPLI 
waveform  and  the  black  lines  enclosing  the  blue  line  represent  the  confidence  interval. 

All  four  pairs  included  a  frontal  and  parietal  electrode  (Fz  and  CP3,  FC3  and  Pz,  F8  and  P02,  and 
C3  and  Pz).  All  waveforms  follow  a  similar  downward  trend  from  3  s  before  the  trigger  pull  to 
the  trigger  pull.  However,  the  first  pair  (2a),  contains  a  significant  point  grouping  immediately 
before  the  trigger  pull,  while  the  other  three  pairs  (2b-d)  contain  significant  groupings  around  3  s 
before  the  trigger  pull.  Another  similarity  is  the  frequency  bands  in  which  the  significant  point 
groupings  occur,  as  three  of  the  four  longest  groups  occur  in  the  theta  (4-8  Hz)  band.  Figure  5c 
is  the  only  significant  point  grouping  that  occurs  in  the  delta  (1-4  Hz)  band.  Further  analyses  are 
needed  to  clarify  the  significance  of  these  pairs,  although  these  particular  pairs  reflect  the  general 
pattern  of  communication  between  parietal  and  frontal  regions  shown  in  figure  4. 
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(a)  (b) 


(c) 


(d) 


Figure  5.  (a)  The  WPLI  analysis  depicted  was  conducted  in  the  theta  band,  4-8  FIz,  on  the  Fz  and  CP3 

electrode.  Thirty  significant  points  were  found  from  -290  to  0  ms  before  the  trigger  pull,  (b)  the 
WPLI  analysis  depicted  was  conducted  in  the  theta  band,  4-8  Hz,  on  the  FC3  and  Pz  electrode. 
Thirty  significant  points  were  found  from  -2990  to  -2700  ms  before  the  trigger  pull,  (c)  the  WPLI 
analysis  depicted  was  conducted  in  the  delta  band,  1-4  Hz,  on  the  F8  and  P02  electrode.  Thirty-two 
significant  points  were  found  from  -2910  to-2600  ms  before  the  trigger  pull,  and  (d)  the  WPLI 
analysis  depicted  was  conducted  in  the  theta  band,  4-8  Hz,  on  the  C3  and  Pz  electrode.  Sixty-four 
significant  points  were  found  from  -2990  to  -2360  ms  before  the  trigger  pull. 

There  are  several  ways  to  improve  on  these  preliminary  results.  Primarily,  all  18  participants  in 
the  study  must  be  analyzed.  The  patterns  observed  in  the  single  participant  strongly  support  the 
extension  of  these  analyses  to  all  participants.  Second,  the  statistical  approach  may  have  been 
suboptimal.  As  was  discussed  regarding  figure  5,  all  four  waveforms  followed  a  generally 
similar  pattern,  yet  the  statistical  approach  revealed  significance  at  different  times  relative  to  the 
trigger  pull.  Alternative  approaches  conducted  in  all  18  participants  may  reveal  a  clearer  result. 
Finally,  this  report  does  not  address  the  potential  issue  of  brain  activity  that  is  correlated  to 
muscle  activity,  which  may  be  the  case  for  the  activity  of  the  F7  electrode  due  to  the  electrode’s 
proximity  to  the  eye.  However,  further  research  is  necessary. 
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2.5  Conclusions 


WPLI  analysis  suggests  that  the  measure  may  be  able  to  capture  task-related  functional  networks 
during  cognitive  tasks,  even  in  the  midst  of  large-scale  movement  artifacts.  Further  research 
looking  across  more  individuals  and  at  more  tasks  is  needed  to  provide  more  definitive  results. 
However,  the  results  from  this  one  participant  support  that  conducting  pair-wise  electrode 
analysis  in  addition  to  individual  electrode  analysis  may  provide  more  information  than 
individual  electrode  analysis  alone.  While  EEG  time-frequency  reveals  a  symmetric  distribution 
of  electrodes  with  significance,  WPLI  analysis  finds  statistical  interdependencies  more  heavily 
lateralized  to  the  left  hemisphere.  In  addition,  WPLI  reveals  communication  with  occipital 
electrodes,  whereas  these  individual  channels  are  not  marked  significant  in  the  time-frequency 
analysis.  The  preliminary  findings  reported  here  suggest  that  the  F7  electrode  may  be  a  hub  of 
communication,  with  several  significant  pairings  with  other  electrodes.  Therefore,  the 
methodology  used  in  this  project  shows  that  WPLI  may  provide  functional  connectivity  data  not 
accessible  in  time-frequency  data.  WPLI  provides  an  interesting  method  for  investigating  how 
the  brain  performs  cognitive  tasks  even  with  large  artifacts  in  the  signal,  allowing  for  further 
investigation  and  understanding  of  Army  tasks  performed  in  complex  environments. 


3.  Project  2:  Functional  Connectivity  of  EEG  Data  Using  the  Weighted 
Phase  Lag  Index  (WPLI):  Robustness  to  Artifacts  in  Real-world 
Environments 


3.1  Authors 

The  authors  on  this  project  are  Scott  Kerick  from  ARL/HRED,  Manuel  Vindiola  from  Dynamics 
Research  Corporation  (DRC),  and  Kaleb  McDowell  from  ARL/HRED. 

3.2  Introduction 

Ecological  psychologists  have  long  known  that  human  behavior  cannot  be  understood 
independent  of  the  complexity  that  the  real-world  imposes  and  that  it  is  critical  to  study  human 
behavior  within  context  (Gibson,  1977;  1979).  As  researchers  have  attempted  to  comprehend  the 
relationship  between  brain  function  and  behavior,  similar  arguments  have  been  made  stressing 
the  importance  of  understanding  brain  function  within  naturalistic  environments  (Gevins,  et  al., 
1995;  Makeig  et  al.,  2009;  Oie  and  McDowell,  2011;  Parasuraman,  2003;  Parasuraman  and 
Rizzo,  2007).  Some  have  gone  as  far  as  proposing  the  idea  that  laboratory-based  studies  isolated 
from  real-world  experience  may  generate  fundamental  misunderstandings  of  the  underlying 
principles  being  investigated  (cf.,  Kingstone  et  al.,  2008).  It  would  thus  seem  important  to 
conduct  neuroscience  investigations  in  natural  settings  to  address  these  concerns.  However,  the 
dominant  brain  imaging  methods  of  the  day  require  participants  to  lie  or  sit  completely  still 
within  large  chambered  machines  (fMRI)  or  under  a  hood  (MEG),  which  imposes  a  substantially 
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different  environment  than  that  experienced  by  humans  performing  tasks  in  natural  settings.  This 
brings  into  question  the  external  validity  and  generalizability  of  much  of  neuroimaging  research 
(Kingstone,  Smilek,  and  Eastwood,  2008;  Kingstone  et  al.,  2003).  Recent  findings  have  bolstered 
these  arguments  by  demonstrating  the  importance  of  context  in  the  differential  engagement  of 
neural  resources  during  task  performance  (Foerde  et  al.,  2006),  and  further  stressed  that  to  better 
understand  how  the  human  brain  operates  in  the  real  world,  considerably  more  research  is 
needed  to  investigate  the  potentially  important  differences  between  the  brain  function  and 
behavior  observed  in  simplistic  and  controlled  laboratory  environments  and  that  observed  in 
complex  and  real-world  settings. 

One  of  the  many  reasons  that  more  neuroscience  has  not  been  conducted  in  real-world  settings  is 
the  challenge  of  recording  quality  signals  from  the  brain  in  natural  environments.  Part  of  this 
challenge  is  due  to  the  characteristics  of  the  technologies  available  for  the  purpose — many  of 
which  are  not  suitable  for  recording  brain  dynamics  in  real-world  environments.  At  present,  the 
only  two  brain  imaging  technologies  that  are  portable  and  wearable  and  thus  capable  of 
recording  brain  dynamics  as  participants  interact  relatively  freely  in  the  naturalistic  environments 
are  EEG  and  functional  near-infrared  spectroscopy  (fNIR).  While  fNIR  is  a  relatively  new 
technology,  EEG  has  been  around  since  the  1920s  (Berger,  1929)  and  over  the  past  five  years, 
there  has  been  a  dramatic  increase  in  commercially  available  and  prototype  mobile  EEG  systems 
(Lin  et  al.,  2010).  Both  technologies  provide  potentially  useful,  yet  different,  indices  of  brain 
activity.  fNIR  provides  a  measure  of  changes  in  blood  oxygen  concentrations  and  blood  volume 
in  response  to  activity  of  neurons  (Izzetoglu  et  al.,  2007),  whereas  EEG  measures  fluctuations  in 
electrical  fields  generated  by  the  net  electrochemical  activity  of  populations  of  pyramidal  cortical 
neurons  and  inter-neurons  (Nunez  and  Srinivasan,  2006).  While  both  of  these  imaging  modalities 
can  provide  relatively  high  quality  signals  in  highly  controlled  settings,  they  both  have  very  low 
signal-to-noise  ratios  in  real-world  settings  due  to  large  artifacts  that  can  be  induced  by 
biological  and  non-biological  sources. 

In  addition  to  the  methodological  challenges  of  obtaining  high  quality  brain  imaging  data  from 
participants  in  natural  environments,  another  major  challenge  is  to  analyze  and  interpret  brain 
function  within  these  environments.  The  brain  is  a  complicated  network  of  neural  elements  that 
has  been  argued  to  function  in  order  to  meet  two  diverse  requirements:  the  local  specialization  of 
neurons  for  processing  specific  types  of  information  and  the  integration  of  information  across  a 
more  global  systems  level  (Friston,  1994;  Rubinov  and  Spoms,  2010;  Spoms,  Tononi,  and 
Edelman,  2000).  While  much  work  has  focused  on  the  local  specialization  of  neural  assemblies, 
it  has  been  argued  that  additional  focus  on  the  systems  level  is  required  in  order  to  advance  our 
current  understanding  of  brain  dynamics  (Bomer,  Sanyal,  and  Vespignani,  2007;  Le  van  Quyen, 
2003;  Varela  et  al.,  2001).  Network  approaches,  and  specifically,  the  inferences  of  functional 
connectivity  upon  which  these  approaches,  in  large  part,  are  built,  and  they  offer  a  potential 
approach  towards  such  a  systems  view.  However,  significant  issues  have  been  identified  for 
deriving  functional  connectivity  from  EEG  data  in  particular.  Of  primary  concern  for  EEG-based 
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methods  are  issues  associated  with  volume  conduction  and  the  use  of  active  reference  electrodes, 
both  of  which  cause  an  artificial  inflation  of  measures  of  association  between  pairs  of  electrodes 
(Bullmore  and  Bassett,  2010;  Stam,  Nolte,  and  Daffertshofer,  2007).  In  an  attempt  to  address  the 
joint  problems  of  volume  conduction  and  common  source  distribution  at  the  scalp,  Stam  et  al. 
(2007)  introduced  the  PLI  as  a  novel  means  of  indexing  the  interdependencies  among 
multichannel  EEG  data.  More  specifically,  the  use  of  the  PLI  is  “. . .  based  upon  the  idea  that  the 
existence  of  a  consistent,  nonzero  phase  lag  between  two  times  series  cannot  be  explained  by 
volume  conduction  from  a  single  strong  source  and,  therefore,  renders  true  interactions  between 
the  underlying  systems  rather  likely”  (Stam  et  al.,  2007,  p.  1179).  Using  simulated  and  actual 
EEG  and  MEG  data,  Stam  et  al.  (2007)  demonstrated  that  PLI  is  less  affected  by  volume 
conduction  and  common  sources  than  are  traditional  measures  of  functional  connectivity  such  as 
phase  coherence  and  the  imaginary  component  of  coherence.  However,  a  limitation  of  PLI  is  the 
discontinuity  of  phase  lag  estimates  between  zero  and  infinitesimally  small  non-zero  phase  lags. 
To  address  this  concern,  Vinck  et  al.  (2011)  introduced  the  WPLI,  which  weights  the  observed 
phase  leads  and  lags  by  the  magnitude  of  the  imaginary  component  of  the  cross-spectrum.  As 
Vinck  et  al.  (201 1)  demonstrated,  WPLI  exhibited  reduced  sensitivity  to  uncorrelated  noise 
sources  and,  concomitantly,  increased  statistical  power  for  detecting  phase  synchronization  as 
compared  with  PLI.  They  further  argue  that  because  WPLI  is  based  on  the  imaginary  part  of  the 
cross-spectrum,  it  is  generally  decreased,  not  increased  with  the  addition  of  noise.  Although  PLI 
and  WPLI  have  been  tested  using  simulated  data  (controlled,  known-truth)  and  on  data  from 
patients  experiencing  epileptic  seizures  (a  condition  exhibiting  drastic  brain  synchronization 
detectable  using  various  methods),  there  exists  a  need  to  validate  these  phase  synchronization 
approaches  on  noisy  data  recorded  in  operational  environments.  Such  validation  could  lead  to  the 
identification  of  a  measure  of  functional  brain  dynamics  that  is  sensitive  to  phase 
synchronization  of  brain  source  signals  but  resilient  to  artifacts  and  would  open  the  door  to 
network  analyses  and  other  systems-level  approaches  to  understanding  brain  function  within 
complex,  naturalistic  environments. 

3.2.1  Objective 

In  response  to  the  challenges  presented  above  and  the  potential  of  WPLI  to  become  a  core 
measure  for  research  in  complex,  artifact-inducing  environments,  we  aim  to  further  understand 
WPLI  with  respect  to  noise  levels  on  the  scale  of  real-world  artifacts.  To  accomplish  this 
objective,  we  examine  EEG  recorded  from  participants  shooting  at  pop-up  targets  in  a  simulated 
environment  with  realistic  weapons  including  requirements  for  dynamic  movements  and  recoil 
immediately  following  weapon  triggering.  From  these  data,  WPLI  was  derived  and  examined  on 
data  epochs  consisting  of  three  types  of  events:  non-event  locked  artifacts  (i.e.,  primarily  non¬ 
brain  activity),  non-event  locked  alpha  bursts  (i.e.,  primarily  brain  activity),  and  trigger-locked 
activity  (i.e.,  a  mix  of  brain  activity  and  artifact  induced  by  recoil).  Similar  to  previous  efforts, 
the  measures  of  functional  connectivity  were  derived  among  all  pair-wise  combinations  of 
electrodes. 
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3.3  Methods  or  Approach 

3.3.1  Participants 

Seventeen  male  U.S.  Marines  and  one  U.S.  Army  Ranger  volunteered  as  participants  (mean  age 
19.8  ±4.2).  All  had  completed  basic  training  and  shooting  qualification  (8  marksmen,  5 
sharpshooters,  5  experts).  All  but  one  of  the  participants  was  right-hand  and  right-eye  dominant; 
one  participant  was  left-hand  and  left-eye  dominant.  The  voluntary,  fully  informed  consent  of  the 
persons  used  in  this  research  was  obtained  as  required  by  32  CFR  219  (1999)  and  AR  70-25 
(1990).  As  is  described  below,  data  for  this  analysis  were  selected  from  12  of  these  original  18 
participants. 

3.3.2  Task  and  Event  Extraction 

The  EEG  data  analyzed  in  this  report  were  acquired  from  participants  during  performance  of 
eight  different  simulated  threat-engagement  shooting  tasks,  an  arithmetic  task,  and  two  resting 
tasks  (for  more  details,  see  Kerick  et  al.,  2007).  For  this  methodological  paper,  we  visually 
inspected  and  manually  selected  three  types  of  transient  events  (200-2000  ms  duration)  from  the 
continuous  EEG  data:  (1)  non-event-related  movement  artifacts  (“Random”),  (2)  spontaneous 
alpha  bursting  (“Alpha  Bursts”),  and  (3)  event-related  brain  activity  with  artifacts  (“Trigger”). 
Random  artifacts  were  primarily  muscle  artifacts  due  to  shifting  positions  between  trials, 
spontaneous  alpha  bursting  events  were  selected  based  on  the  onset  of  clearly  visible  alpha 
oscillatory  patterns  observable  across  multiple  channels,  and  event-related  brain  activity  with 
artifacts  were  selected  based  on  existing  trigger  event  markers  in  the  data,  which  were  always 
immediately  followed  by  weapon  recoil  and  mechanical/movement  artifacts.  Figure  6  shows  an 
example  event  of  each  type  within  continuous  EEG  recordings.  Fifty  trials  of  each  event  type 
were  extracted  by  epoching  the  data  around  each  event  from  -1  s  to  2  s,  where  0  s  corresponded 
to  the  onset  of  artifact  contamination  for  Random  events,  the  onset  of  alpha  oscillations  for 
Alpha  Bursting  events  and  the  time  of  the  trigger  pull  for  the  Trigger  events.  The  150  total 
selected  events  were  sampled  from  12  of  the  participants  and  10  of  the  conditions. 
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1.  Non-Event-Related  Artifact 
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Figure  6.  Example  illustrating  three  types  of  events  extracted  from  continuous  EEG 
channel  data.  (Top  to  bottom):  (1)  non-event-related  artifact,  (2)  transient 
alpha  bursting,  and  (3)  event-related  artifact.  Note  that  these  signals  are  from 
data  prior  to  epoching  and  that  a  longer  than  3  s  time  window  is  illustrated. 


3.3.3  EEG  Acquisition 

EEG  data  were  acquired  at  500  Hz  using  a  Neuroscan  system  (NuAmps  40-channel)  with  34 
active  electrodes,  2  horizontal  and  2  vertical  eye  channels,  and  2  mastoid  reference  channels  (left 
and  right).  The  right  mastoid  served  as  a  common  reference  for  all  other  channels  during  online 
recording.  Online  filter  settings  were  DC- 100  Hz. 

3.3.4  EEG  Pre-processing 

All  pre-processing  was  executed  using  EEGLAB  (ver.  9.0.2.2b;  Delorme  and  Makeig,  2004) 
running  on  Matlab  (2008a)  installed  on  a  64-bit  Linux  operating  system.  Raw  continuous  EEG 
data  from  each  condition  were  bandpass  filtered  from  1-50  Hz  using  a  two-way  least-squares 
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finite  impulse  response  (FIR)  filter  (using  EEGLAB  function  eegfilt.m).  After  bandpass  filtering 
the  data,  files  for  each  subject  and  condition  were  saved  for  use  in  two  separate  data  processing 
streams.  For  the  first  processing  stream,  data  were  concatenated  from  each  of  the  1 1  conditions 
for  each  subject  separately  and  then  epoched  into  non-overlapping  500-ms  windows.  Bad 
channels  were  identified  and  removed  from  the  filtered,  concatenated,  epoched  data,  by  visual 
inspection  aided  by  automatic  channel-rejection  procedures  based  on  kurtosis  and  joint 
probability  functions  (using  EEGLAB  function pop_rejchan.m).  Once  bad  channels  were 
removed,  automatic  artifact-detection  procedures  were  applied  (using  EEGLAB  function 
pop_rejmenu.m )  based  on  thresholds  for  kurtosis,  joint  probability,  abnormal  distribution,  and 
abnormal  trends,  as  well  as  visual  inspection  to  remove  any  epochs  consisting  of  non¬ 
stereotypical  transient  artifacts  (e.g.,  gross  muscle  contamination,  cable  sway  artifacts,  etc.).  The 
filtered,  artifact-reduced,  concatenated,  epoched  data  (i.e.,  “clean  EEG”)  were  then  re-referenced 
to  an  average  reference  derived  from  all  active  channels  (excluding  eye  channels).  The  rationale 
for  concatenating  data  across  all  conditions  and  removing  non-stereotypical  artifacts  was  to 
provide  the  cleanest  and  greatest  amount  of  data  available  for  ICA  decomposition  (Onton  and 
Makeig,  2006).  Extended  infomax  ICA  (Bell  and  Sejnowski,  1995;  Makeig,  Bell,  Jung,  and 
Sejnowski,  1996)  was  then  applied  to  decompose  the  clean  EEG  data  into  spatially  fixed, 
maximally  temporally  independent  components  (using  EEGLAB  function  runica.m )  and  dipole 
models  were  fit  to  each  of  the  independent  components  (ICs;  using  EEGLAB  plug-in  DipFit2.2 ). 
ICs  were  then  classified  into  three  categories  based  on  dipole  source  models,  topographic  maps, 
spectra  and  activation  waveforms:  (1)  clearly  brain-source,  (2)  clearly  artifact-source  (these  IC 
reflected  primarily  eye  blinks  and  muscle  artifacts),  and  (3)  mixed-source.  However,  while  the 
ICs  derived  using  these  methods  can  be  used  to  separate  brain-source  from  artifact-source  under 
good  conditions,  large  transient  artifacts  permeate  nearly  all  ICs.  The  outputs  from  the  ICA 
decomposition  (IC  weights,  inverse  weights,  sphering  matrix,  component  indexes  and  IC 
classifications)  and  dipole  source  models  were  then  saved. 

For  the  second  processing  stream,  we  retrieved,  for  each  participant,  each  of  the  band-pass- 
filtered  (1-50  Hz)  data  files  from  each  of  the  1 1  conditions  and  removed  the  same  bad  channels 
removed  in  the  first  processing  stream  and  then  re-referenced  to  the  average  of  all  remaining 
good  channels  (however,  the  previously  identified  bad  epochs  were  not  removed).  The  outputs 
from  the  ICA  decomposition  and  dipole  source  models  were  then  copied  into  the  data  structures 
for  each  participant  and  each  of  the  1 1  conditions.  ICs  classified  as  artifact-source  were  removed 
and  the  EEG  data  were  then  back-projected  to  channels  void  of  those  ICs  resembling  artifacts 
(Jung  et  al.,  2000;  Onton  et  al.,  2006)  (i.e.,  “ICA-filtered  EEG  data”).  Note  that  we  did  not 
remove  any  of  the  transient,  non-stereotypical  artifacts  in  any  of  the  1 1  condition  data  sets, 
beyond  those  that  were  removed  by  ICA  decomposition  (e.g.,  eye-blinks  [ICs  1,  2,  5,  10]  and 
some  muscle  artifacts  [ICs  21  and  22];  see  figure  7)  as  the  purpose  of  this  study  was  to  observe 
the  effects  of  such  artifacts  on  WPLI. 
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Figure  7.  Channel  power  by  event  and  frequency  band  for  each  of  the  50  individual  trials.  Within  each  plots, 
trials  are  sorted  by  the  length  of  the  EEG  event  (shorter  events  on  top).  Three  of  the  six  frequency 
bands  (theta,  alpha  II,  gamma)  are  plotted  in  columns.  All  three  event  types  (Random,  Alpha,  and 
Trigger)  are  plotted  in  rows.  A  black  bar  below  each  plot  shows  when  the  distributions  of  power 
values  were  significantly  higher  than  baseline  (-1  Os). 


3.3.5  Weighted  Phase  Lag  Index  (WPLI) 

WPLI  was  derived  from  the  continuous  ICA-filtered  EEG  channel  data  for  each  participant  in 
each  condition.  WPLI  was  derived  for  each  pair-wise  combination  of  channel  data.  Below  we 
describe  how  WPLI  was  derived.  Data  were  first  band-pass  filtered  using  a  4th  order  Butterworth 
filter  at  the  following  frequency  bands:  delta  (1-4  Hz),  theta  (4-8  Hz),  alpha  I  (8-10  Hz),  alpha 
II  (10-13  Hz),  beta  (13-30  Hz),  and  gamma  (30-45  Hz).  The  Hilbert  transform  was  then  applied 
to  each  of  the  band-pass  filtered  data  to  obtain  an  analytical  signal.  WPLI  is  defined  as  follows: 


wPLiJi^m 

E{\Z{X}\} 


(1) 


where  3  [X  J  is  the  imaginary  part  of  the  cross-spectrum  (Vinck  et  al.,  201 1).  WPLI  estimates 
were  computed  over  four  cycles  for  each  frequency  band  (i.e.,  varying  length  windows  were 
defined  as  four  times  the  length  of  one  cycle  of  the  lower  bound  of  the  band-pass  filter  and  were 
stepped  in  10  ms  increments  [five  samples]  for  each  frequency  band). 


24 


For  each  signal  type  (spectral  band  power,  WPLI),  event  type  (Random,  Alpha  Burst,  Trigger), 
and  band-pass  fdter  setting  (delta  through  gamma),  separate  significance  tests  were  conducted 
using  the  Wilcoxon  rank  sum  test  for  equal  across-trial  distributions  between  pre-event  data  and 
each  time  step  through  the  entire  length  of  the  trial.  That  is,  for  each  of  the  36  possible  time 
series,  301  (time  points  from  [-1  2])  rank  sum  tests  were  performed  comparing  distributions 
between  the  pre-event  data  (100  time  points  by  50  trials)  and  each  specific  time  point  (1  time 
point  by  50  trials).  Significance  was  determined  using  Bonferroni  corrections  (p  <  2.3  x  106)  to 
account  for  the  chance  occurrence  that  any  of  the  several  tests  were  significant.  In  addition,  at 
least  five  consecutive  time  points  had  to  be  significant  in  order  to  be  included  in  the  results.  Note 
that  PLI  measures  were  also  computed  for  each  bandpass  filter  setting  and  for  each  condition.  As 
nearly  identical  results  were  observed  between  WPLI  and  PLI,  only  WPLI  will  be  reported. 

3.4  Results  and  Discussion 

The  challenge  of  recording  quality  brain  signals  in  complex,  real-world  environments  is  a 
significant  obstacle  to  neuroscience  and  related  fields.  Here,  we  examine  a  novel  measure  of 
functional  brain  connectivity  (i.e.,  WPLI  [Vinck  et  al.,  2011])  that  shows  promise  for  being  a 
robust  measure  of  brain  communication  in  the  face  of  real-world,  non-brain  artifacts.  In  order  to 
assess  the  robustness  of  WPLI  to  realistic,  real-world  artifact  contamination  of  EEG  recordings, 
three  types  of  events  representing  non-brain  artifact,  brain  activity,  and  a  mix  of  non-brain 
artifact  and  brain  activity  were  examined.  To  confirm  the  severity  of  the  artifact,  advanced  signal 
processing  of  the  EEG  data  (i.e.,  ICA;  Makeig,  Bell,  Jung,  and  Sejnowski,  1996)  was  performed 
that,  under  typical  laboratory  conditions,  can  isolate  brain-sources  from  non-brain  sources.  Both 
non-brain  artifact  (Random)  and  mixed  (Trigger)  events  were  selected  for  which  the  artifact 
overpowered  the  source  discrimination.  We  argue  here  that  the  results  support  that  WPLI  is 
insensitive  to  large  non-brain  artifacts  and  appears  to  be  sensitive  to  brain  activity. 

3.4.1  Spectral  Power 

Figure  7  illustrates  the  spectral  power  for  individual  trials  in  the  channel  data,  respectively,  for 
three  frequency  bands.  The  black  bars  below  the  individual-trial  plots  in  figure  7  as  well  as  the 
information  in  table  1  indicate  the  time  points  that  are  significantly  higher  than  baseline  [-1  0  s]. 


Table  1.  Latency  ranges  (low-high  ms)  of  significant  power  increase  for  each  event 
type  and  frequency  band. 


Delta 

Theta 

Alpha  I 

Alpha  II 

Beta 

Gamma 

Random 

-10:160 

180:260 

-80:300 

-80:320 

-90:300 

-20:400 

-30-450 

Alpha 

- 

- 

-30:1590 

-60:1460 

160:350 

- 

Trigger 

-140:330 

-120:400 

-140:680 

-100:290 

-60:310 

-40:340 
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Random 


The  channel  data  exhibited  significant  (p  <  2.3  x  106)  power  increases  relative  to  baseline  in  each 
of  the  six  frequency  bands.  These  significant  increases  occurred  from  90  to  10  ms  prior  to  onset 
and  extended  for  260  to  450  ms  after  event  onset  (see  table  1).  The  data  indicate  a  high  level  of 
broadband  power  (e.g.,  compare  scales  across  all  three  event  types),  which  is  consistent  with 
artifact.  While  significant  power  increases  across  all  six  frequency  bands  existed,  on  a  trial-by- 
trial  basis  power  increases  following  event  onset  appear  to  be  less  variable  in  the  higher 
frequency  bands,  which  is  also  consistent  with  muscle-generated  artifact  that  generally  has  more 
power  at  the  higher  end  of  the  EEG  spectrum. 

Alpha  Bursting 

The  channel  data  exhibited  significant  (p  <  2.3  x  106)  power  increases  relative  to  baseline  in  the 
two  alpha  bands  and  the  beta  band.  The  significant  increases  generally  occurred  from  -60  to 
-20  ms  prior  to  onset  and  extended  for  350  to  1590  ms  after  event  onset  (see  table  1).  Overall, 
the  data  indicate  higher  pre-  and  post-event  power  in  lower  frequency  bands  (e.g.,  compare 
apparently  random  power  increases  in  theta  to  very  low  power  in  gamma),  which  is  consistent 
with  cortically-generated  activity.  On  a  trial-by-trial  basis,  the  power  distribution  across 
frequency  bands  appears  to  have  similar  variability.  These  results  demonstrate  the  properties  of 
neurally-generated  signal  without  substantial  artifact. 

Trigger 

The  data  show  significant  (p  <  2.3  x  106)  power  increases  relative  to  baseline  in  each  of  the  six 
frequency  bands.  These  significant  increases  occurred  from  -180  to  -40  ms  prior  to  onset  and 
extended  for  330  to  500  ms  after  event  onset  (see  table  1).  In  this  event,  brain  activity  had  to 
occur  prior  to  the  trigger  pull.  The  relatively  early  onsets  of  the  power  increases  occurred  well 
before  the  recoil  (i.e.,  large  movement  artifacts),  and  this  timeframe  supports  the  interpretation 
that  the  brain  activity  was,  in  part,  contributing  the  power  in  this  event.  The  data  also  indicate 
similar  broadband  power  levels  as  those  observed  for  the  alpha  bursting  trials  (e.g.,  compare 
scales)  and  indicate  similar  temporal  patterns  of  power  distribution  across  frequency  bands.  We 
interpret  the  relatively  similar  levels  of  power  across  frequencies  to  reflect  a  combination  of  both 
the  lower  frequency  dominance  that  is  typical  of  brain  activity  combined  with  the  higher 
frequency  dominance  that  is  typical  of  muscle-generated  artifact.  On  a  trial-by-trial  basis,  the 
post-event  power  increases  around  trigger  onset  appear  very  consistent  across  frequency  bands, 
which  may  reflect  the  consistency  in  both  brain  activity  and  recoil  artifacts  associated  with  a 
trigger  pull. 

3.4.2  Weighted  Phase  Lag  Index 

In  discussing  the  WPLI  results,  we  need  to  consider  several  points.  First,  WPLI  is  a  phase-based 
measure  that  generally  will  only  decrease  when  artifact  is  evident  in  the  signal  (Vinck  et  al., 
2011).  Second,  here  we  examined  averaged  WPLIs,  which  should  enable  the  ability  to  detect  a 
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signal  decrease  associated  with  large  artifacts;  however,  at  the  same  time  it  will  limit  the  ability 
to  uncover  the  characteristics  of  complex  brain-generated  network  activity  (which  we  leave  for 
future  efforts).  Third,  the  windowing  to  compute  the  WPLI  smears  the  signal  across  time,  which 
makes  it  critical  to  not  only  consider  the  timing  of  WPLI  effects  relative  to  brain/non-brain 
activity,  but  also  if  the  effects  are  asymmetrical  to  that  activity. 


Figure  8  illustrates  WPLI  for  individual  trials  for  three  frequency  bands,  respectively.  The  red 
bars  below  the  individual-trial  plots  in  figure  8  as  well  as  the  information  in  table  2  indicate  the 
time  points  for  which  WPLI  is  significantly  different  than  baseline  [-1  0  s].  Black  bars  are  also 
included  below  each  individual-trial  plot  to  reference  significant  time  periods  observed  for  the 
spectral  power. 
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Figure  8.  WPLI  by  event  and  frequency  band  for  each  of  the  50  individual  trials.  Within  each  plot,  trials  are  sorted 
by  the  length  of  the  EEG  event  (shorter  events  on  top).  Three  of  the  six  frequency  bands  (theta,  alpha  II, 
gamma)  are  plotted  in  columns.  All  three  event  types  (Random,  Alpha,  and  Trigger)  are  plotted  in  rows.  A 
black  bar  below  each  plot  shows  when  the  distributions  of  power  values  were  significantly  higher  than 
baseline  [-10  s].  A  red  bar  below  each  plot  shows  when  the  distributions  of  WPLI  were  significantly 
higher  than  baseline  [-1  0  s]. 
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Table  2.  Latency  ranges  (low-high  ms)  of  significant  WPLI  increase  for  each  event  type  and 
frequency  band. 


Delta 

Theta 

Alpha  I 

Alpha  II 

Beta 

Gamma 

Random 

- 

- 

- 

- 

- 

10:110 

Alpha 

- 

- 

-50:960 

1300-1650 

-60:1040 

- 

- 

Trigger 

- 

30:780 

180:680 

240:410 

-90:340 

-40:20 

80:330 

Random 

The  EEG  data  revealed  a  significant  (p  <  2.3  x  106)  increase  in  WPLI  in  the  gamma  band  only. 
The  significant  WPLI  increase  began  10  ms  after  event  onset  and  extended  for  1 10  ms  (see 
table  2).  Both  in  the  time  period  of  significance  and  on  a  trial-by-trial  basis,  WPLI  does  not 
appear  to  reflect  the  artifact-related  power  increases  (e.g.,  compare  duration  of  significance  in  the 
EEG  channel  data  gamma  band  for  power  [black  bar]  and  WPLI  [red  bar]).  Further,  the  short 
periods  of  WPLI  change  are  increasing  not  decreasing  as  would  be  expected  due  to  noise  artifact. 
It  is  conceivable  that  the  WPLI  changes  are  actually  reflecting  some  level  of  neural  activity  that 
could  be  associated  with  muscle  activity,  even  though  this  activity  was  not  “event-locked”  to  a 
cognitive  activity  (e.g.,  see  trigger  pull).  Importantly,  however,  the  results  here  illustrate  very 
small  changes,  especially  when  considering  the  amplitude  of  the  artifact  and  the  strength  of  the 
results  in  the  following  two  event  types. 

Alpha  Bursts 

The  EEG  channel  data  revealed  a  significant  (p  <  2.3  x  106)  increase  in  WPLI  in  the  alpha  I  and 
II  bands  only.  In  the  EEG  channel  data,  the  significant  increases  occurred  60  ms  (alpha  I)  and 
50  ms  (alpha  II)  prior  to  event  onset  and  extended  for  940  ms  (alpha  I)  and  1040  ms  (alpha  II) 
after  event  onset;  there  was  also  a  later  period  of  significance  in  the  alpha  I  band  (see  table  2). 
These  data  support  that  WPLI  is  capable  of  reflecting  brain-based  functional  connectivity. 
Further,  both  in  the  time  periods  of  significance  and  on  a  trial-by-trial  basis,  the  WPLI  values  of 
the  channel  data  do  not  precisely  coincide  with  the  brain-related  alpha  power  increases  (e.g., 
compare  duration  of  significance  in  the  EEG  channel  data  alpha  II  band  for  power  [black  bar] 
and  WPLI  [red  bar]).  While  supporting  that  WPLI  reflects  functional  connectivity  in  this  case, 
these  data  also  suggest  a  level  of  independence  between  amplitude  generation  and  synchronized 
global  activity  which  has  been  previously  supported  (Makeig  2002;  2003). 

Trigger 

The  EEG  channel  data  revealed  a  significant  (p  <  2.3  x  106)  decrease  in  WPLI  in  all  frequency 
bands  except  delta.  The  significant  decreases  generally  occurred  90  to  40  ms  before  trigger  pull 
in  the  higher  frequency  bands  (beta  and  gamma)  and  30  to  240  ms  after  the  trigger  pull  in  middle 
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frequency  bands  (theta,  alpha  I  and  alpha  II)  and  extended  for  330  to  780  ms  after  the  trigger 
pull.  These  results,  particularly  taken  together  with  that  of  the  Random  and  Alpha  Burst  events, 
suggest  that  the  WPLI  significances  are  not  reflecting  muscle-generated  artifact:  the  level  of 
significance  for  the  trigger  events  are  substantially  longer  and  across  broader  frequencies  than 
that  of  the  random  artifact  events  even  though  the  artifacts  were  much  smaller,  and  there  exists 
large  asymmetrical  discrepancies  in  time  period  between  the  channel  WPLI  and  the  channel 
spectral  data.  Again  while  supporting  that  WPLI  reflects  functional  connectivity  in  the  Trigger 
events,  the  data  also  suggest  some  independence  between  amplitude  generation  and  synchronized 
phase  activity  (Makeig  2002;  2003). 

3.5  Conclusions 

In  summary,  we  have  presented  a  computational  approach  that  addresses  some  of  the  many 
issues  and  challenges  associated  with  acquiring  and  interpreting  brain  imaging  data  in  the  form 
of  EEG  signals  recorded  in  an  operationally-relevant,  non-traditional  laboratory  environment. 
Namely,  we  applied  a  novel  measure  of  functional  brain  connectivity  based  on  phase 
synchronization  (i.e.,  WPLI  [Vinck  et  al.,  2011])  among  multiple  recording  sites  to  investigate 
the  performance  of  this  measure  in  the  face  of  volume-conducted-artifacts  in  EEG  signals,  which 
present  a  pervasive  problem  for  translational  neuroscience  research  and  understanding  brain 
function  in  real-world  environments.  Together,  the  identification  of  a  measure  of  functional  brain 
dynamics  that  is  sensitive  to  phase  synchronization  of  brain  source  signals  but  resilient  to 
artifacts  would  potentially  provide  important  new  insights  into  how  the  brain  functions  in 
complex,  real-world  settings.  That  is,  the  novel  WPLI  measure  appears  simultaneously  largely 
insensitive  to  the  large  artifacts  found  in  realistic,  real-world  environments  and  sensitive  to  brain¬ 
generated  functional  connectivity.  These  properties  support  that  WPLI  has  a  strong  potential  to 
become  not  only  a  core  measure  for  research  into  understanding  brain  function  in  complex, 
artifact-inducing  environments,  but  may  also  become  a  core  component  in  the  development  of 
brain-computer  interaction  technologies.  However,  before  the  potential  of  WPLI  can  be  fully 
realized,  a  substantial  amount  of  future  research  is  necessary  to  uncover  the  specific  sensitivities 
of  WPLI  to  brain  dynamics. 


4.  Project  3:  Validating  Functional  Connectivity  Measures  Using  Neural 
Mass  Models 


4.1  Authors 

The  authors  on  this  project  are  Manuel  Vindiola  from  DRC,  Stephen  Gordon  from  DCS 
Corporation,  Jean  Vettel  from  ARL/HRED,  and  Kaleb  McDowell  from  ARL/HRED. 
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4.2  Introduction 

The  human  brain  can  be  conceived  of  as  a  massively  parallel,  massively  distributed  set  of 
electrochemical  processors  called  neurons  interconnected  with  one  another  by  axons. 
Collectively,  trillions  of  neurons  and  billions  of  axons  form  a  complex  information  processing 
network.  The  combined  activity  and  communication  of  these  neurons  underlie  human  cognition 
and  behavior. 

When  observed  individually,  each  neuron  fires  according  to  a  non-stationary,  random  Poisson 
point  process;  however,  when  a  population  of  neurons  is  combined,  their  average  firing  rates 
form  a  continuous,  periodic,  oscillatory  pattern  of  neuronal  activity.  Oscillatory  neural  activity 
has  been  identified  in  numerous  brain  functions  including  attention,  navigation,  perceptual 
binding,  and  memory  (Gray  et  al.,  1989).  Oscillatory  brain  rhythms  can  be  phase-coupled  across 
distant  brain  regions.  Two  regions  are  phase-coupled  or  phase-synchronized  when  they  display  a 
systematic  phase-offset.  Phase  synchronization  has  been  identified  as  a  way  to  create  a  flexible 
communication  structure  between  distant  brain  regions  (Engel  et  al.,  2005;  Fries,  2005). 

Communication  among  brain  regions  form  networks  of  brain  activity,  creating  stable  brain  states 
that  capture  a  global  level  of  neural  processing  thought  to  subserve  ongoing  task  demands  and 
behavioral  performance.  Task-relevant  brain  networks  dynamically  reconfigure,  capturing  the 
temporal  dynamics  of  the  task.  These  ongoing  network  reconfigurations  may  engage  different 
brain  regions,  alter  the  nature  of  the  communication  between  regions,  etc.  In  order  to  measure 
and  understand  the  nature  of  these  ongoing  brain  dynamics,  neuroimaging  analysis  methods  have 
been  developed  to  capture  communication  between  pairs  of  recorded  brain  signals. 

EEG  is  a  method  of  measuring  electrical  activity  on  the  scalp,  and  evidence  suggests  that  EEG 
signals  preserve  the  oscillatory  characteristics  of  the  firing  patterns  in  the  neural  regions  that 
generate  them.  Computational  and  statistical  advancements  over  the  last  decade  have  led  to  the 
development  of  a  number  of  functional  connectivity  measures  that  can  be  applied  to  EEG 
signals.  There  are  several  classes  of  connectivity  measures,  including  ones  based  on  correlation, 
auto-regression,  and  phase  synchronization.  In  this  analysis,  we  focus  on  three  undirected,  phase- 
based  functional  connectivity  measures  that  have  been  designed  to  identify  pair-wise  phase 
synchronizations  among  separate  EEG  signals:  imaginary  component  of  Coherence  (ImC)  (Nolte 
et  al.,  2004),  debiased  Weighted  Phase  Locking  Index  (dWPLI)  (Vinck  et  al.,  2011),  and  phase¬ 
locking  value  (PLV)  (Lachaux  et  al.,  1999).  These  three  methods  were  chosen  because  they  offer 
three  different  views  into  an  issue  that  typically  confounds  connectivity  analysis  in  EEG  data — 
volume  conduction. 

Volume  conduction  has  two  major  components:  (1)  brain  sources  that  are  mixed  together  as  they 
pass  through  the  scalp,  and  (2)  artifact  sources  on  the  scalp  (e.g.,  eye  movements,  muscle 
contractions).  In  typical  analyses  of  brain  function,  the  aim  is  to  separate  artifact  from  brain  and 
unmix  the  brain  sources.  These  three  measures  show  different  sensitivities  to  these  sources  of 
volume  conduction.  In  general,  PLV  is  negatively  affected  by  both  brain  and  artifact  sources  of 
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volume  conduction,  while  both  ImC  and  dWPLI  minimize  the  effect  of  artifact  sources  by 
requiring  that  the  two  signals  have  a  systematic  phase  relationship  that  occurs  at  a  temporal 
delay.  If  the  two  channels  are  instantaneously  in  phase,  then  ImC  and  dWPLI  both  return  low 
functional  connectivity  values,  under  the  premise  that  removing  zero-phase  volume  conducted 
noise  signals  improves  the  functional  connectivity  analysis  even  if  some  real  communications 
between  brain  regions  are  sacrificed.  ImC  will  also  return  a  low  connectivity  value  for  phase- 
lagged  synchronizations  if  the  brain  signals  and  non-brain  artifacts  induce  large  amplitudes  in  the 
signal.  Unlike  ImC  and  dWPLI,  PLY  returns  high  functionally  connectivity  values  so  long  as  the 
two  signals  are  phase-synchronized  and  the  delay  is  consistent.  Since  we  know  volume 
conduction  introduces  a  linear  mixing  of  the  signals,  i.e.,  instantaneous  phase  synchronization,  it 
is  worth  examining  whether  the  methods  that  attempt  to  deal  with  this  effect  perform  equally 
well,  or  better,  in  the  more  ideal,  unmixed  case.  Are  these  functional  connectivity  measures  able 
to  detect  task-relevant  changes  in  ongoing  brain  network  dynamics?  To  address  this  question,  we 
applied  these  functional  connectivity  measures  to  simulated  oscillatory  signals  generated  by  a 
pair  of  coupled  regions  in  a  neural  mass  model  in  order  to  evaluate  how  well  the  measures 
detected  changes  in  ongoing  brain  network  dynamics  when  volume  conduction  was  not  present. 
If  the  measures  can  detect  task-relevant  brain  dynamics,  they  provide  a  critical  capability  to 
measure  and  monitor  brain  function  in  Army-relevant  tasks  and  settings  and  thus  open  an  avenue 
for  neurotechnology  development  to  improve  Soldier-system  performance. 

4.2.1  Objective 

In  this  project,  a  simulation  was  developed  to  validate  the  functional  connectivity  measures;  that 
is,  can  the  measure  recover  the  network  dynamics  when  temporal  dynamics  of  phase 
synchronizations  are  known  to  exist?  Here,  we  have  a  neural  mass  model  with  several  regions. 
Each  region  is  parameterized  to  oscillate  with  a  specific  frequency  distribution  and  is  coupled 
with  other  regions  to  form  a  network.  We  begin  with  the  simplest  case  of  two  regions  with  each 
region  oscillating  with  one  of  two  possible  frequency  distributions  connected  to  each  other  with 
one  of  two  possible  network  configurations.  We  first  confirm  that  the  simulated  regions  generate 
the  oscillatory  signals  expected  from  a  given  network,  and  then  we  run  each  of  the  three  phase 
connectivity  measures  to  evaluate  what  network  dynamics  are  captured. 

4.3  Methods  or  Approach 

4.3.1  Neural  Mass  Models 

The  neural  mass  models  (NMMs)  were  constructed  according  to  David  and  Friston’s  (2003) 
description.  Each  region  within  the  NMM  is  designed  to  reproduce  a  time  signal  with  amplitude 
and  oscillation  frequencies  similar  to  those  observed  in  local  field  potential  (LFP)  recordings 
obtained  by  measuring  the  combined  post  synaptic  potential  of  hundreds  to  thousands  of  real 
neurons  in  a  cortical  column.  A  key  component  of  these  NMMs  that  makes  them  suitable  for  this 
study  is  that  both  the  signal  patterns  they  produce  and  those  observed  in  EEG  recordings  are 
oscillatory.  A  sample  signal  from  a  region  in  a  NMM  is  presented  in  figure  9. 
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Figure  9.  The  time  course  of  a  single  region  in  a  NMM  in  isolation  with  two 
internal  subpopulations.  The  output  is  generated  by  nonlinearly 
combining  a  weighted  mixture  of  80%  from  a  subpopulation  oscillating 
at  1 1  FIz  and  20%  from  a  subpopulation  oscillating  at  43  Hz. 

A  NMM  region  is  composed  of  a  number  of  internal  subpopulations,  where  each  subpopulation 
is  parameterized  to  oscillate  at  one  specific  frequency.  The  output  signal  of  a  neural  mass  region 
is  constructed  by  computing  a  nonlinear  sigmoid  transformation  over  a  linear  weighted 
combination  of  each  of  the  subpopulations.  The  subpopulations  create  oscillatory  patterns 
through  the  interaction  of  excitatory  and  inhibitory  components  of  the  subpopulation.  Each 
subpopulation,  which  is  loosely  based  on  the  structural  organization  of  cortical  macro-columns, 
contains  three  layers:  a  pyramidal  layer,  an  excitatory  spiny  stellate  layer,  and  an  inhibitory 
intemeuron  layer.  The  cortical  output  is  modeled  by  a  summed  weighting  of  all  the  excitatory 
pyramidal  layers.  The  pyramidal  layer  receives  input  from  the  excitatory  spiny  stellate  and 
inhibitory  intemeuron  layers  of  the  region.  The  excitatory  spiny  stellate  layer  receives  input  from 
the  pyramidal  layers  of  other  neural  mass  regions  that  are  part  of  the  same  functional  network 
(neighbors),  external  subcortical  bottom-up  inputs  modeled  by  a  stochastic  Gaussian  process, 
and  feedback  from  the  pyramidal  layers  of  the  same  region.  The  inhibitory  intemeuron  layer 
receives  excitatory  feedback  input  from  the  pyramidal  layers  of  the  same  region,  and  it  sends 
inhibitory  output  to  these  same  pyramidal  layers.  Figure  10  depicts  a  schematic  representation  of 
the  connectivity  among  the  three  parts  of  a  subpopulation  and  the  connectivity  between  two 
different  neural  mass  regions. 
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Figure  10.  Architecture  of  two  interconnected  neural  mass  regions.  Each  region  is  composed  from  three 

columnar  layers.  The  excitatory  spiny  stellate  layer  receives  the  external  input  coming  from  other 
regions.  The  pyramidal  layer  sends  output  to  other  regions.  The  interaction  between  the  excitatory 
spiny  stellate  and  inhibitory  interneuron  layers  give  rise  to  oscillatory  activation  patterns  in  the 
pyramidal  layer.  Each  region  may  have  multiple  subpopulations,  each  of  which  drives  oscillations  at  a 
different  fundamental  frequency.  In  this  figure,  region  1  has  three  internal  subpopulations  and  region 
2  has  two.  O  is  the  weighted  sum  of  the  sub-populations  and  S  is  the  sigmoid  function. 

4.3.3. 1  Neural  Mass  Region  State  Space 

Formally,  the  NMM  el  is  a  dynamical  system  composed  of  six  coupled  first-order  differential 
equations  per  subpopulation  in  each  region.  These  equations  were  solved  numerically  using  the 
fourth-order  Runge-Kutta  method  from  the  SPM8  (Wellcome  Department  of  Cognitive 
Neurology,  University  College  London,  London,  UK)  toolbox  for  Matlab  (2011a,  The 
Mathworks.  Natick,  MA.). 

The  equations  for  the  state  space  are 

Excitatory  Spiny  Stellate 


v  i 

u,sp,  1 


=  X 


u,sp,  1 
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where  u  is  a  neural  mass  region,  sp  is  a  subpopulation  in  region  u,  N  is  the  number  of 
subpopulations  in  region  u,  He  tunes  the  maximum  excitatory  post  synaptic  potential,  xe  models 
properties  of  the  dendritic  tree,  Networku  is  the  input  to  region  u  from  other  regions  in  the 
functional  network,  External  is  the  input  from  sub-cortical  areas,  Co  is  the  average  number  of 
synaptic  contacts  ,  and  S  is  a  sigmoid  function  that  transform  average  membrane  potential  into  an 
average  rate  of  action  potentials  fired  by  the  neurons  in  the  region.  S  takes  the  form 
2c 

S(inp )  = -  0  ,  where  eo=2.5,  r  =  .56,  vo  =  6.  wu  /(  is  the  input  weighting  for  subpopulation  k, 

1  +e^vo  _inp^ 

pu  is  the  average  input  from  subcortical  units,  and  pu (t)  is  the  deviations  from  the  distribution 
N(  pu  ,22).  kw  ll  is  the  coupling  weight  for  input  from  region  w  to  region  u.  S(  yH.(?  -  S))  is  the 
standard  deviation  of  the  pyramidal  output  from  region  w  at  time  t  —  S.  For  these  simulations, 

5  =  10  ms  and  all  standard  deviations  were  computed  between  [t-500,  t-S].  kwu  is  a  coupling 

parameter  that  is  recomputed  at  each  time  step  to  ensure  that  the  variance  of  the  combined  input 
from  Internal  and  External  is  conserved. 
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H;  tunes  the  maximum  inhibitory  post  synaptic  potential,  i;  models  properties  of  the  dendritic 
tree;  all  other  parameters  represent  the  same  properties  already  described. 

Pyramidal  Cells 
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Average  Firing  Rate  for  Region 
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4.3.2  Simulation  Structure 

The  networks  were  constructed  by  coupling  two  neural  mass  regions  together.  Each  region  of  the 
pair  was  instantiated  with  an  1 1-  and  43 -Hz  internal  subpopulation.  The  weighting  of  the 
subpopulations  was  either  (1)  80%  of  the  1 1-Hz  subpopulation  and  20%  of  43-Hz  subpopulation, 
which  yielded  a  broadly  tuned  oscillator  with  a  peak  in  the  power  spectrum  at  10  Hz,  or  (2)  20% 
of  the  1 1-Hz  subpopulation  and  80%  of  the  43-Hz  subpopulation,  which  yielded  a  broadly  tuned 
oscillator  with  a  peak  in  the  power  spectrum  at  40  Hz.  All  possible  combinations  of  10-  and 
40-Hz  units  were  connected  using  unidirectional  and  bidirectional  couplings  for  a  total  of  seven 
pair  combinations.  Figure  11  depicts  all  simulated  pair  combinations  that  were  instantiated. 


Ten  simulated  trials  were  performed  for  each  pair  combination.  Each  trial  lasted  for  a  total 
duration  of  1 1  s  with  a  sample  rate  of  1  ms.  The  first  second  was  included  to  allow  state-space 
stabilization,  and  it  was  discarded  prior  to  computing  the  power  spectrum  and  functional 
connectivity  measures.  The  connection  between  the  units  was  active  between  1  s  through  3  s  and 
5  s  through  8  s.  The  connection  was  inactive  at  all  other  times. 

4.3.3  Functional  Connectivity  Measures 

Each  of  the  functional  connectivity  measures  were  computed  as  shown  in  the  following 
equations: 

IMC  =  ^=EI{Z]Z^  (12) 

^{|Z12|}£{|Z22|} 
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(14) 


dWPLI  = 


i  E&v&m 

E{\Z{ZX}\) 


where  Z/  and  Z2  are  the  Fourier  spectra  of  the  two  signals  being  compared,  (*)  denotes  complex 
conjugate,  and  61  and  62  are  the  relative  phases  for  each  signal  determined  by  representing  Z  as 
the  phasor  Z  =  Ael6.  It  should  be  noted  that  the  expectation  operator  E{. }  is  taken  over  trials,  thus 
the  results  shown  were  obtained  by  averaging  over  ten  separate  trials  for  each  simulated 
condition. 

To  compute  the  functional  connectivity  measures  over  time,  a  moving  time  window  of  500  ms 
was  used  with  an  incremental  step  of  25  ms.  This  window  size  was  chosen  empirically  to  provide 
optimal  tradeoff  between  time  and  frequency  resolution.  Within  each  time  window,  Fourier 
spectra  Z/  and  Z2  were  computed  using  the  Matlab’s  fast  Fourier  transform  method  resulting  in 
2000  equally  sized  frequency  bands  spanning  a  frequency  range  of  0  to  1000  Hz.  The  value  2000 
was  determined  by  setting  the  nonequispaced  fast  Fourier  transform  (NFFT)  parameter  in  Matlab 
to  be  four  times  the  sample  rate  (1000  Hz),  which  meant  the  original  signal  had  to  be  zero- 
padded  to  accurately  compute  the  desired  number  of  frequency  components.  For  presentation 
purposes,  the  results  have  been  limited  to  a  range  of  1  to  100  Hz  as  this  contained  all  the 
oscillatory  components  of  interest. 

4.4  Results  and  Discussion 

4.4.1  Result  1:  Validation  of  the  Neural  Mass  Models 

The  first  set  of  results  seeks  to  establish  that  the  NMMs  produce  oscillatory  firing  patterns 
appropriate  for  use  with  phase-based  functional  connectivity  measures.  In  order  to  establish  this, 
two  properties  of  the  generated  signals  were  investigated.  The  signal  time  course  was  examined 
for  oscillatory  patterns  sensitive  to  whether  the  coupling  between  regions  was  on  or  off.  When 
the  coupling  is  on,  the  amplitude  of  the  receiving  region  should  show  a  large  increase  over  when 
the  coupling  is  off.  This  basic  result  should  be  present  because  the  receiving  region  receives  a 
larger  input  when  the  coupling  is  on  and  this  increase  should  be  reflected  in  the  output. 

The  other  expected  property  is  a  difference  in  the  frequency  power  spectrum  distribution  (PSD) 
when  the  coupling  is  on  and  off.  The  exact  nature  of  the  change  depends  on  the  properties  of  the 
two  regions  and  the  type  of  connection  that  couples  them.  A  unidirectional  coupling  should  show 
a  change  in  the  PSD  of  the  receiving  region  and  an  unchanged  PSD  in  the  sending  region.  When 
the  coupling  is  in  both  directions,  both  regions  should  show  a  change  in  the  PSD. 

This  section  presents  results  from  four  of  the  seven  coupled  pairs  that  were  simulated.  The  first 
two  coupled  pairs  are  identical  regions  coupled  in  both  directions:  one  pair  is  a  10-Hz  region 
coupled  to  another  10-Hz  region,  and  the  other  pair  is  a  40-Hz  region  coupled  to  another  40-Hz 


region.  The  time  course  activity  plotted  in  figure  12  shows  that  both  regions  oscillate  with  an 
average  depolarization  between  7.4  and  7.6  when  the  coupling  is  off.  When  the  coupling  is 
turned  on  between  1  s  through  3  s  and  5s  through  8s,  they  both  show  increases  in  the  average 
depolarization.  The  10  Hz  <— >  10  Hz  pair  oscillated  with  an  average  depolarization  between  7.2 
and  7.9,  and  the  40  Hz  40  Hz  pair  oscillated  between  6.8  and  8.89.  These  changes  indicate 
that  the  NMMs  are  sensitive  to  the  communication  between  regions  and  show  a  measurable 
change  during  coupling. 


Figure  12.  Time  course  and  PSD  plots  for  two  pairs  of  identical  bicoupled  neural  mass  regions.  Time  courses  show 
that  they  produce  oscillatory  firing  patterns,  and  the  firing  rate  increases  greatly  when  the  coupling  is  on 
(1-3  s  and  5-8  s).  Power  spectral  density  plots  show  that  both  units  have  the  same  power  distribution 
both  when  the  coupling  is  off  and  on.  When  coupling  is  on,  both  regions  show  a  small  downshift  in  the 
peak  frequency  as  well  as  a  large  increase  in  power  at  that  frequency  and  its  second  harmonic. 

Figure  12  shows  that  the  PSD  plots  for  both  regions  of  a  pair  are  virtually  identical  when  the 
coupling  is  off;  both  units  of  the  10  Hz  10  Hz  pair  have  identical  unimodal  PSDs  with  a 
peak  at  10  Hz  and  a  1  to  50  Hz  distribution.  Both  regions  of  the  pair  share  identical  bimodal 
PSDs  when  the  coupling  is  on,  but  the  shape  of  the  distribution  differs  from  the  distribution 
when  coupling  is  off.  The  central  peak  shifts  down  a  small  amount  to  9.4  Hz  when  the  coupling 
is  on  and  has  a  substantial  power  increase  compared  to  when  the  coupling  is  off.  A  second  power 
increase  is  also  present  at  a  frequency  close  to  the  second  harmonic  of  central  peak. 

The  40  Hz  4--^  40  Hz  pair  shows  a  similar  pattern  to  the  10  Hz  4--^  10  Hz  pair.  The  two 
regions  have  identical  unimodal  PSDs  centered  at  40  Hz  when  the  coupling  is  off.  When  the 
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coupling  is  turned  on,  the  shape  of  the  PSD  changes  and  this  change  occurs  in  both  of  the 
regions.  Compared  to  the  PSD  when  the  coupling  is  off,  the  new  PSD  shows  a  small  downshift 
in  the  peak  frequency  as  well  as  a  substantial  increase  in  spectral  power  at  that  frequency  and  its 
second  harmonic.  Both  the  10  Hz  < — >  10  Hz  and  40  Hz  <-->  40  Hz  pairs  also  show  a  large 
broadband  spectral  power  increase  with  a  1-Hz  peak.  The  reason  for  this  low  frequency  power 
increase  is  still  under  investigation,  but  it  may  be  due  to  induced  synchronizations  that  can  be 
seen  in  the  amplitudes  during  coupling. 

The  second  two  coupled  pairs  are  a  40-Hz  region  coupled  in  only  one  direction  to  another  40-Hz 
region,  and  10-  and  40-Hz  regions  coupled  to  each  other  in  both  directions.  When  the  coupling  is 
off,  the  time  series  and  PSD  of  the  40-  and  10-Hz  regions  are  the  same  as  those  in  figure  12. 
When  the  coupling  is  on,  the  sending  40-Hz  region  in  the  40  Hz->40  Hz  pair  (red  lines  in  the  top 
row  of  figure  13)  does  not  show  any  difference  from  when  the  coupling  is  off  in  either  its  time 
series  or  in  the  PSD  plots.  Conversely,  the  receiving  region  (blue  lines  in  the  top  row  of 
figure  13)  shows  a  similar  pattern  to  the  40-Hz  units  in  the  40  Hz  ^ — >  40  Hz  pair.  There  are 
significant  amplitude  increases  in  the  time  series  and  a  slight  down  shift  and  large  power 
increase  of  the  peak  frequency,  its  second  harmonic,  and  at  1  Hz.  Neither  the  amplitudes  nor  the 
powers  are  as  strong  as  they  were  in  the  bidirectionally  coupled  variant.  There  are  no  changes 
between  this  region  and  two  regions  in  40  Hz  4-  40  Hz  pair  at  the  timeframes  when  the 

coupling  is  off;  this  is  expected  since  this  region  is  not  receiving  any  input. 


Figure  13.  Time  course  and  power  spectral  density  plots  for  two  pairs  of  NMMs.  Time  courses  show  that 
they  produce  oscillatory  firing  patterns  and  the  firing  rate  increases  greatly  when  the  coupling  is 
on  (1-3  s  and  5-8  s).  Power  spectral  density  plots  show  that  the  power  distribution  is  different 
when  coupling  is  off  compared  to  when  it  is  on. 
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When  the  coupling  is  on,  the  10  Hz  < — >  40  Hz  pair  produces  the  now  expected  pattern  of 
results.  Both  units  oscillated  with  the  same  PSD  observed  in  the  other  network  variants  when  the 
coupling  was  off  and  had  bimodal  PSDs  with  peaks  at  1  Hz,  the  expected  frequency,  and  its 
second  harmonic  when  coupling  was  on  (figure  13,  bottom  row).  The  other  expected  result, 
which  can  be  seen  clearly  in  figure  14,  is  that  both  regions  show  a  slight  downward  shift  in  the 
frequencies.  The  10-Hz  region’s  PSD  peak  shifted  from  10  to  8  Hz  and  the  40-Hz  region’s  PSD 
peak  shifted  from  40  to  35  Hz.  The  reason  for  the  shift  is  still  under  investigation,  but  this  result 
suggests  that  a  pair  of  coupled  regions  can  cause  each  other  to  oscillate  at  different  frequencies 
than  they  would  in  isolation.  Results  of  this  type  provide  a  potentially  valuable  insight  for  results 
obtained  in  EEG  recordings.  A  neural  region  might  intrinsically  oscillate  at  a  different  frequency 
than  is  observed  when  it  is  processing  as  part  of  a  functional  network.  It  is  also  possible  that  the 
frequency  at  which  a  given  region  oscillates  changes  as  a  function  of  its  functional  network.  This 
result  supports  such  a  possibility  because  it  provides  preliminary  evidence  that  the  topography  of 
the  functional  network  exerts  an  influence  on  a  region’s  oscillation  frequencies.  Future  work  is 
necessary  to  examine  the  extent  to  which  this  manifests  itself  in  larger  networks  as  well  as 
whether  specific  network  motifs  affect  the  likelihood  of  a  frequency  shift. 


Power  Spectral  Density  Power  Spectral  Density 
Coupling  Off  Coupling  On 


Figure  14.  Power-frequency  plots  for  10  Hz  40  Hz  pair  when  coupling  is  off  and  on.  When  the  coupling  is  off, 
both  units  have  a  unimodal  power  distribution  centered  at  their  respective  dominant  frequencies.  When 
coupling  is  on  both  units  have  a  bimodal  power  distribution.  One  peak  comes  from  the  unit’s  intrinsic 
oscillation  rate  and  the  other  arises  because  influences  from  the  other  unit.  The  most  interesting  part  is 
that  the  coupling  causes  a  frequency  downshift  in  both  units.  The  10-Hz  unit  has  peaks  at  8  and  35  Hz 
and  the  40-Hz  unit  has  power  peaks  at  3,  7,  35,  and  66  Hz. 


A  new,  but  predicted,  result  when  coupling  is  on  is  that  the  PSDs  of  both  regions  contain  a 
combination  of  each  region’s  intrinsic  frequency  as  well  as  an  additional  component  coming 
from  the  other  region.  For  example,  the  10-Hz  unit’s  PSD  shows  a  large  spectral  power 
component  at  ~10  Hz  and  a  second  large  component  at  ~35  Hz.  This  new  component,  which  is 
not  present  when  the  coupling  is  off,  reflects  input  from  the  40-Hz  unit.  The  same  is  true  of  the 
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PSD  of  the  40-Hz  unit.  In  addition  to  its  intrinsic  frequency,  it  also  contains  an  8-Hz  component 
coming  from  the  10-Hz  unit  when  the  coupling  is  on  between  the  pair.  Both  of  these  results 
contain  the  expected  power  spectrum  components  from  communication  between  regions  with 
different  intrinsic  frequencies. 

Collectively,  these  results  confirm  that  the  NMM  simulation  produces  activations  that  are 
affected  by  the  oscillation  frequencies  and  connectivity  of  two  region  networks.  The  simulated 
data  thus  provides  an  avenue  to  validate  functional  connectivity  measures  and  examine  if  the 
measures  recover  the  time  course  of  paired  communication  within  the  network.  The  simulated 
functional  activity  from  each  of  these  networks  was  analyzed  with  three  undirected,  phase-based 
functional  connectivity  measures  (ImC,  dWPLI,  and  PLV).  Two  interesting  results  from  the 
measure  comparison  are  discussed  in  section  4.4.2. 

4.4.2  Result  2:  Limitation  of  ImC  and  dWPLI 

One  situation  in  which  PLV  provides  a  better  estimate  of  functional  connectivity  than  ImC  or 
dWPLI  is  when  two  identically  parameterized  regions  are  bidirectionally  coupled.  With  this 
configuration,  neither  ImC  nor  dWPLI  will  identify  the  phase  synchronizations  induced  by  the 
functional  connectivity.  Figure  15  demonstrates  that  PLV  is  high  when  the  coupling  is  on, 
whereas  ImC  and  dWPLI  return  consistently  low  connectivity  values  regardless  of  whether 
coupling  is  on  or  off. 


PIV  ImC  dWPLI 


Figure  15.  Phase-based  functional  connectivity  estimates  using  PLV,  ImC  and  dWPLI.  Only  PLV  is  able  to  identify 
30-50  Hz  phase  synchronization  but  the  measure  also  has  a  strong  tendency  to  incorrectly  identify  phase 
synchronizations  at  lower  frequencies  regardless  of  whether  the  coupling  is  on  or  off.  ImC  and  dWPLI 
fail  to  identify  phase  synchronizations  at  any  frequency. 
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When  two  identical  oscillators  are  coupled  together,  they  very  quickly  synchronize  their 
oscillations  in-phase  or  anti-phase.  Our  results  indicate  that  this  is  true  of  the  NMMs  even  though 
there  is  a  small  time  delay  of  10  ms  between  the  two  regions.  ImC  and  dWPLI  return  low  values 
with  this  configuration  because  they  were  specifically  constructed  to  do  this  when  two  signals 
oscillate  in-phase  or  anti-phase.  Both  were  designed  to  be  resistant  to  noise  artifacts  from  volume 
conduction,  where  a  zero-phase  relationship  is  discarded  as  artifact;  however,  as  figure  15  shows, 
they  also  miss  legitimate  in-phase  synchronizations  produced  by  interacting  brain  regions. 
Interestingly  if  the  magnitude  of  the  value  is  ignored  and  its  stability  across  time  is  examined 
instead,  there  is  a  clear  and  consistent  value  throughout  the  entire  duration  of  communication. 
PLV,  which  does  not  minimize  zero-phase  oscillations,  does  identify  the  in-phase 
synchronization  that  occurs  when  the  coupling  is  on,  but  it  is  expected  that  PLV  will  be  more 
adversely  affected  by  motion  and  environmental  artifacts  when  tested  on  empirical  EEG  datasets 
collected  with  real-world  movements  and  noise. 

4.4.4  Result  3:  dWPLI  Has  Greatest  Sensitivity  in  Low  Frequency  Bands 

All  three  measures  successfully  identify  phase  synchronizations  in  the  dominate  20-60  Hz 
frequency  range  for  the  40  Hz  40  Hz  pair  and  the  25 — 40  Hz  frequency  range  for  the 
10  Hz  40  Hz  coupled  pair  presented  in  figure  16.  dWPLI  is  more  sensitive  with  phase 
synchronization  in  the  lower  frequency  bands  than  PLV  and  ImC.  The  40->40  coupled  pair  in 
figure  16  (top  pair)  shows  that  dWPLI  identifies  a  low  frequency  phase  synchronization  when 
the  coupling  transitions  from  off  to  on  and  again  when  it  transitions  back  off.  PLV  also  identifies 
low  frequency  phase  synchronizations,  but  the  problem  is  that  PLV  identifies  so  many  false 
synchronizations  when  the  coupling  is  off  that  is  it  is  difficult,  if  not  impossible,  to  identify  the 
true  low  frequency  phase  synchronizations  at  coupling  transitions.  ImC,  on  the  other  hand,  is  not 
sensitive  to  the  lower  frequency  synchronization  at  all;  this  is  evidenced  by  the  fact  that  the  value 
is  near  zero  across  the  entire  time  duration.  Another  curiosity  for  this  pair  is  that  the  ImC  and 
dWPLI  values  for  this  pair  are  consistently  zero  at  40  Hz.  This  is  because  of  the  same  reason  the 
value  is  zero  for  bidirectionally  coupled  pairs;  during  communication  these  two  regions  are 
oscillating  in-phase  or  anti-phase  with  each  other. 
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Figure  16.  (D)  Phase-based  functional  connectivity  estimates  using  PLV,  ImC  and  dWPLI.  All  three  measures 

identify  20-60  Hz  phase  synchronizations  when  coupling  is  on  but  in  the  lower  1-5  Hz  frequency  range, 
only  dWPLI  identifies  a  reliable  difference  in  phase  synchronization. 

The  10  Hz  40  Hz  pair  (figure  16,  bottom)  provides  additional  support  for  this  claim.  In  this 
case,  dWPLI  identifies  sustained  low  frequency  phase  synchronizations  that  are  not  reliably 
identified  by  PLV  or  ImC.  Once  again  the  phase  synchronizations  identified  by  PLV  cannot  be 
trusted  because  it  also  identifies  many  low  frequency  phase  synchronizations  at  points  when  the 
coupling  is  off.  ImC’s  performance  is  slightly  better  than  with  the  previous  pair.  It  picks  up  on 
some  phase  synchronization  in  the  3  Hz  band  when  the  coupling  is  on  but  it  is  still  much  less 
sensitive  than  dWPLI.  In  summary,  all  three  measures  can  capture  the  30-60  Hz  phase 
synchronizations,  but  only  dWPLI  reliably  identifies  1-5  Hz  phase  synchronization. 

4.5  Conclusions 

This  project  implemented  oscillatory  NMMs  in  order  to  simulate  oscillatory  functional  activity, 
and  then  simulated  data  from  these  models  were  used  to  validate  three  undirected,  phase-based 
functional  connectivity  measures  (ImC,  dWPLI,  and  PLV). 

The  first  set  of  results  established  that  the  NMMs  produce  oscillatory  signals,  and  when  coupled, 
their  output  is  influenced  by  the  input  coming  from  paired  neighbors  during  fixed  windows  of 
communication.  That  is,  the  NMMs  successfully  produced  oscillatory  behavior  and  directed 
communication  within  the  network.  Consequently,  the  model  appears  to  be  a  good  tool  to 
evaluate  if  functional  connectivity  measures  can  capture  the  temporal  dynamics  of  network 
communication. 

The  second  set  of  results  indicated  that  all  three  connectivity  measures  capture  the  temporal 
dynamics  of  network  communication  for  higher  frequencies  (i.e.,  30-60  Hz)  when  nodes  in  the 
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network  oscillate  at  different  frequencies.  The  measures  revealed  different  sensitivities,  however, 
when  the  network  communication  occurred  at  the  same  frequency  or  the  communication 
occurred  at  lower  frequencies  (1-5  Hz).  First,  only  PLV  could  capture  communication  between 
bicoupled  nodes  oscillating  at  the  same  frequency,  since  both  ImC  and  dWPLI  were  designed  to 
minimize  contribution  of  zero-phase  signal  synchronizations.  This  indicates  the  tradeoff  between 
resilience  to  volume  conduction  noise  in  the  signal  and  sensitivity  to  one  type  of  oscillatory 
communication  between  network  regions.  Second,  dWPLI  captured  communication  at  the  lower 
frequencies  within  the  networks  most  robustly.  Here,  the  value  of  the  PLV  measure  was  too 
noisy  in  the  lower  frequencies  to  distinguish  between  coupled  and  uncoupled  timeframes,  while 
ImC  was  rather  insensitive  at  this  lower  range.  Combined,  these  results  have  begun  to  identify 
the  relative  strengths  and  weaknesses  of  each  of  the  measures,  and  it  suggests  that  this  modeling 
effort  holds  promise  to  investigate  the  sensitivities  of  the  connectivity  measures  in  controlled, 
simulated  networks  to  better  understand  how  to  interpret  results. 

The  second  set  of  results  also  indicated  that  the  topography  of  the  network  had  a  significant 
impact  on  the  performance  of  the  functional  connectivity  measures.  In  the  10  Hz  4--^  40  Hz 
pair,  the  uncoupled  PSD  revealed  the  expected  peaks  for  the  intrinsic  frequencies  (10  and 
40  Hz,  respectively),  but  during  the  time  windows  where  the  two  nodes  are  coupled,  the  PSDs 
for  the  nodes  changed,  shifting  the  10  to  8  Hz  and  the  40  to  35  Hz.  In  other  words,  the  properties 
of  the  functional  network  modulated  the  node’s  intrinsic  oscillation  frequency,  and  this  could 
have  important  consequences  for  interpretation  of  the  measure  given  known  associations  in  the 
literature  between  particular  frequency  ranges  and  cognitive  constructs.  The  hope  is  that 
developing  an  understanding  of  how  functional  connectivity  measures  are  affected  by  network 
structure  in  a  simulation  where  the  underlying  connectivity  is  known  will  provide  guidance  and 
insight  for  selecting  which  measures  to  use  on  real  human  data  when  the  functional  and 
structural  networks  are  unknown. 

Critically,  future  efforts  will  expand  the  model  to  simulate  larger  networks  with  more  diverse 
topologies  and  topographies  in  order  to  better  link  the  simulation  work  with  empirical  research. 
The  functional  networks  that  have  been  explored  thus  far  are  two  unit  networks  with 
unidirectional  or  bidirectional  connections  between  them.  The  functional  networks  observed  in 
human  and  primate  experiments  have  revealed  functional  networks  consisting  of  dozens  of 
interconnected  regions.  Systematically  examining  network  topology,  interconnectivity,  and 
topology  has  potential  for  helping  with  the  interpretation  of  functional  connectivity  results 
obtained  from  more  complex  networks.  Larger  networks  may  provide  more  unforeseen  insights 
into  the  effect  of  network  configuration  on  functional  connectivity  measures.  Larger  networks 
will  also  allow  for  a  more  thorough  and  systematic  examination  of  how  different  network 
configurations  change  the  intrinsic  oscillations  of  the  regions. 

A  second  important  avenue  for  future  research  is  a  thorough  exploration  of  the  parameter  space. 
Specifically,  it  is  important  to  understand  which  parameters  are  stable  with  respect  to  variation 
and  which  parameters  cause  instabilities  when  changed.  Expanding  to  larger  networks  may  cause 
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the  models  to  become  unstable  when  using  the  current  parameter  values.  Furthermore,  growing 
the  networks  to  dozens  (or  hundreds)  of  regions  may  require  introducing  strong  boundary 
conditions  over  some  parameters.  Systematically  studying  the  sensitivity  of  parameters  to 
variation  will  provide  insight  into  which  parameters  may  require  these  strong  boundary 
conditions.  In  addition,  the  stability  of  a  parameter  may  be  modulated  by  the  size  of  the  network, 
since  some  parameters  that  do  not  play  a  significant  role  in  small  networks  may  become  more 
sensitive  as  the  networks  become  larger. 

In  summary,  the  results  here  indicate  that  the  NMM  simulation  can  reproduce  expected 
oscillatory  behavior  in  two-node  networks,  and  even  these  simple  networks  capture  substantial 
differences  among  the  three  tested  functional  connectivity  measures.  Future  efforts  will  explore 
larger  networks,  parameters  to  keep  those  larger  network  stable,  and  additional  types  and  classes 
of  functional  connectivity  measures.  This  simulation  work  is  expected  to  augment  empirical 
efforts  to  develop  and  validate  functional  connectivity  measures  in  order  to  understand  network- 
based  brain  dynamics  in  complex,  Army-relevant  settings. 


5.  Project  4:  Design  and  Implementation  of  a  Numerical  Technique  to 
Inform  Anisotropic  Hyperelastic  Finite  Element  Models  Using  Diffusion- 
weighted  Imaging 


5.1  Authors 

The  authors  on  this  project  are  Reuben  H.  Kraft  ARL/WMRD  and  Amy  M.  Dagro  ARL/WMRD 
(ARL-TR-5796,  in  press). 

5.2  Introduction 

Anisotropy  is  prevalent  in  numerous  biological  tissues,  which  are  comprised  of  fibers,  or  bundles 
of  cells,  aligned  in  a  uniform  direction.  For  example,  within  the  white  matter  regions  of  the 
human  brain,  axons  tend  to  form  complex  fiber  tracts  that  promote  anatomical  connectivity  and 
communication.  The  bundles  of  axonal  fibers  can  be  visualized  using  noninvasive  tools  such  as 
diffusion  tensor  magnetic  resonance  medical  imaging  (DTI).  The  basic  principle  behind  DTI  is 
that  water  diffuses  more  rapidly  in  the  direction  aligned  with  the  internal  structure  and  more 
slowly  as  it  moves  perpendicular  to  the  preferred  direction.  From  the  diffusion  tensor  computed 
from  the  imaging,  diffusion  anisotropy  measures  such  as  the  FA  can  further  be  computed.  In 
addition,  the  principal  direction  of  the  diffusion  tensor  can  be  used  to  estimate  the  white  matter 
connectivity  of  the  brain.  Visualization  of  the  white  matter  connectivity  is  called  tractography. 

An  example  of  DTI  tractography  is  shown  in  figure  17.  Note  that  DTI  is  one  particular  method  or 
application  of  the  broader  DWI. 
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Figure  17.  Using  DTI,  axonal  fiber  tractography  can  be  overlaid  on  a  finite  element  mesh  and  used  to  construct 
a  transversely  isotropic  model.  Tractography  consists  of  complex  fibers  shown  from  (b)  dorsal, 

(c)  right  lateral  side,  and  (d)  posterior  views. 

The  anisotropic  properties  of  biological  tissue  are  important  to  capture  within  a  computational 
framework  because  the  axonal  fiber  tracts  have  been  reported  to  be  approximately  three  times 
stiffer  than  the  surrounding  matrix  (Arbogast  et  al.,  1999).  In  fact,  numerous  efforts  have  focused 
on  incorporating  details  of  fiber  tractography  into  a  computational  framework  where  many  use 
transversely  isotropic  viscoelastic  constitutive  models  to  represent  the  bulk  and  fiber  mixture. 

For  example,  Ning  et  al.  (2006)  conduct  experiments  and  optimization  techniques  to  develop  a 
transversely  isotropic  viscoelastic  constitutive  model  for  a  porcine  brain  stem  tissue  sample. 
Wright  et  al.  present  an  analytical  and  computational  model  that  also  treats  the  white  matter  as  an 
anisotropic,  hyperelastic  material  using  DTI  to  incorporate  the  structural  orientation  of  the  neural 
axons  into  the  model.  Wright  et  al.  show  that  the  degree  of  injury  that  is  predicted  in  a 
computational  model  of  DTI  is  highly  dependent  on  the  incorporation  of  the  axonal  orientation 
information  and  the  inclusion  of  anisotropy  into  the  constitutive  model  for  white  matter. 

The  work  presented  here  extends  previous  efforts  by  using  the  transversely  isotropic  hyperelastic 
model  formulation,  and  describes  a  new  design  and  implementation  for  informing  the  finite 
element  model  using  DTI  in  a  three-dimensional,  parallel  computing  framework.  Besides 
biological  materials,  the  algorithms  described  in  this  report  can  be  applied  to  a  broad  range  of 
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engineering  materials,  as  well.  For  example,  Kevlar  armor  protection  systems  contain  fiber 
directions  that  rotate  with  the  curvature  of  the  surface.  The  finite  element  implementation 
included  here  sets  the  appropriate  framework  for  incorporating  fiber  directions  into  finite 
element  models  used  in  a  diverse  range  of  applications. 

In  this  report,  we  begin  by  describing  the  continuum  theory  behind  the  implementation  of  the 
transversely  isotropic  hyperelastic  model  in  section  5.3.1  followed  by  results  of  a  validation 
study  of  the  implementation  within  the  finite  element  framework  in  section  5.3.2.  Then,  section 

5.3.3  presents  an  algorithm  for  linking  DTI  and  the  transverse  isotropic  model.  Section  5.3.4 
discusses  an  extension  of  the  research  to  modeling  the  white  matter  tissue  of  the  human  brain  in 
a  three-dimensional  finite  element  code.  Then  section  5.3.5  offers  some  concluding  remarks  and 
opportunities  for  future  work. 

5.3  Methods,  Results,  and  Discussion 


5.3.1  Continuum  Theory  of  Transversely  Isotropic  Hyperelasticity 


5. 3. 1.1  Kinematics  of  Finite  Deformations 


We  begin  by  defining  the  deformation  gradient  defined  as 


F  =  J“1/3F 


(15) 


where  x  is  the  deformed  configuration  of  a  material  particle  and  X  is  the  reference  position.  The 
ratio  of  deformed  to  undeformed  volumes  is  given  by  the  determinant  of  F: 

J  =  det(F)  (le 

Some  materials  behave  differently  in  bulk  and  shear,  so  it  is  most  beneficial  to  perform  a 
multiplicative  decomposition  of  F  into  volume-changing  (dilational)  and  volume-preserving 
(distortional)  parts  (equation  18).  To  do  this,  a  deviatoric  deformation  gradient  is  defined  in 
which  the  volume  change  is  eliminated: 


It  then  follows  that: 


F  =  J”1/3  F 

(17) 

c  =  ftf 

(18) 

B  =  FFT 

(19) 

where  C  and  B  are  referred  to  as  the  modified  right  Cauchy-Green  tensor,  and  the  modified  left 
Cauchy-Green  tensor,  respectively.  The  modified  principle  invariants  of  the  Cauchy-Green 
deformation  tensor  are 
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I 1  -  tr  C  =  tr  B 

h  =  i[(lrC)2-tr(c2)]=i[(trB)2-lr(B2)] 

1 3  ~  det  C  —  det  B 

Two  additional  invariants  are  introduced  in  order  to  capture  the  mechanics  of  transverse  isotropic 
materials: 

74  -  ao-Cao  (23) 

h  =  ao.C2ao  (24) 

where  ?4  and  7s  arise  directly  from  the  anisotropy  and  describe  the  properties  of  the  fiber  family 
and  its  interaction  with  the  other  material  constituents.  The  vector  ao  is  introduced  to  describe  the 
undeformed  fiber  direction  where  we  have  assumed  that  the  only  anisotropic  property  of  the 
solid  comes  from  the  presence  of  the  fibers.  For  a  material  that  is  reinforced  by  only  one  family 
of  fibers,  the  stress  at  a  material  point  depends  not  only  on  the  deformation  gradient  F,  but  also 
on  that  single  preferred  direction,  which  is  called  the  fiber  direction. 

5. 3. 1.2  Decoupled  Representation  of  Strain  Energy 

Following  the  approach  taken  by  others  (Ning  et  al.,  2006;  Holzapfel  et  al.,  2000;  Weiss  et  al., 
1996)  in  the  representation  of  transversely  isotropic  hyperelasticity,  a  unique  decoupled 
representation  of  the  strain  energy  function  per  unit  volume  in  terms  of  the  five  invariants  is 
given  by 


(20) 

(21) 

(22) 


W(/i,  l2y  ^5)  =  ^3)  "h  'ffaniso(74j  7o) 

where  J2J3)  describes  the  response  of  the  isotropic  matrix,  and  describes  the 

directional  contribution  of  the  fiber  reinforcement.  Note  that  for  an  incompressible  material, 

J3  =  .p  =  1 


5. 3. 1.3  Uncoupled  Stress  Tensor 

For  a  hyperelastic  material,  the  2nd  Piola-Kirchhoff  stress  is  derived  from  the  strain  energy  as 


S 


2dC 

2j-2/3DEV 


(C) 

dC 


-pJC~l 


(26) 

(27) 


where  the  operator  DEV  H  -  H  ~  a  (H  :C)C  ’  extracts  the  deviatoric  part  of  a  tensor  in  the  reference 

configuration  and  dJ  is  the  hydrostatic  pressure.  Using  the  chain  rule,  equation  27  can  be 

further  written  as 
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s  =  ‘2J~2/3  DEV 


-pJC 


-1 


dh. 


where  the  terms  oc  are  given  as: 

dh 

ec 

dh 

dC 

du 

dC 

dh 

dC 


=  I 


.0=1  -.W3  '  ^  - 

(28) 

I 

(29) 

7,1  -C 

(30) 

a0  ®ao 

(31) 

^o^C’ao-l-ao’C^ao 

(32) 

where  I  is  a  2nd  rank  identity  tensor. 

A  push-forward  operation  on  the  second  Piola-Kirchoff  stress  tensor  S  to  the  current 
configuration,  and  a  scaling  with  the  inverse  of  the  volume  ratio,  transforms  equation  30  to  the 
Cauchy  stress  tensor  of  a  transversely  isotropic,  compressible  hyperelastic  material: 


a  —  —dev 

u 

T  =r2  ^ 


f2Eft 

8C 


-pi 


^  r  _  _  _ _  _  _  _ 

a  =  -dev  |^('kjl  +  /i'k2)  B  -  ¥2B  +  /4ty4a  <8>  a  +  /4^5  (a  <8>  Ba  +  aB  <g>  a)  -pi 


where 

<M-|  =  H-  '(!•]:  i)  i 

Vi  =  — 

dh 

a  —  Fao 

in  which  a  describes  the  configuration  of  ao  during  the  deformation. 


(33) 

(34) 


(35) 

(36) 

(37) 


5.3. 1.3  Uncoupled  Stress  Tensor  for  a  Particular  'F 

Within  the  framework  of  incompressible,  transversely  isotropic  hyperelasticity,  the  strain  energy 
function  can  be  written  in  the  uncoupled  form  as 


'k(C)  -  Vvoi(J)  +  (c) 


(38) 
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where  'Fvoi(J)  represents  the  volumetric  strain  energy  contribution  and  ^isoc(C)  represents  the 
isochoric  strain  energy  function.  Following  Weiss  et  al.  (1996),  the  isochoric  part  of  the  strain 
energy  can  be  further  subdivided  into  separate  functions  associated  with  the  modified  invariants: 


*(C)  =  'Fvoi(J)  +  Fi  (h,h)  +  F2  (h) 


(39) 


with: 


F  i  =  <?,(/!  — 3)+C2(/2-3) 
F2  =  C-i  [exp  (I4  -  1)  -  74] 


(40) 

(41) 


where  Ci,  C2,  and  C3  are  constants  obtained  from  parameter  fits  to  experimental  data.  Equation 
40  corresponds  to  a  modified  Mooney-Rivlin  model,  and  equation  27  describes  an  exponential 
behavior  in  the  fiber  direction  that  is  a  characteristic  of  most  soft  tissues. 

Then,  similar  to  equation  27,  the  second  Piola-Kirchoff  stress  is  given  by 


S  =  2«7-2/3DEV 


^isoeCC) 

dC 


a*vol(j) 


oj 


JC 


-1 


(42) 


where  the  partial  derivative  of  the  isochoric  part  is  given  by 


^iSoc(C) 


ac 


dVisoc  (C)  ,  d*isoc  (c)  7 1 T  a*isoc  (C)  _  ,  a*isoc  (c) 

- Wr' ' — ”  +  ——=——1 1  1  - - —  + - '= - 


dli  9h  J  dl'2 

so  then  substituting  equation  43  into  the  DEV  operator  in  equation  42 
DEV 


dh 


ao  ®  a0 


(43) 


isoc  (C) 

ac 

aow  (c)  a*  isoc  (C)7 

H - - i  1 


dh 


dl2 


fr£W(C)_  ,  dV  isoc  (C) 

dh  c+ 


a#iS0C{c)_ 

— = — 1 1  + 1 — = — 1 2  -i - — — 1 4 

Oh  dh  dh 


dh 

C  1 


ao  ®  ao 


(44) 


Then,  the  second  Piola-Kirchoff  stress  is  given  by 
S  =  2  J~2/3 


a^soc  (c)  a*  isoc  (C)7\, 
Wi  &h  y I_ 


wtep)=  (c) 

^  + 


dh 


_i  / WteRj  2^oc (c) + a^isoc (c) 

3  \  dli  1  a/2  2  diA  y 


dh 

&KAJ) 

dJ 


a«  ®  ao 
JC"1 


(45) 
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where 


0*vol(</} 

dJ 

=  p 

(46) 

d*iSoc(C) 

=  Ci 

(47) 

0/i 

0*isoc  (C) 

=  a 

(48) 

dh 

S^iSOC  (C) 

-  C3  [exp  (/4  -  l)  -  l] 

(49) 

Then,  substituting  into  equation  45 
S  -  -pJC~l+ 

2 J~2/s  [(Ci  +  C2/i)  I  -  C'i C  +  C3  [exp  (i4  -  l)  -  l]  a«  ®  a,, 
-1/3  (Cji  +  2C272  +  C3  [exp  (l4  -  1)  -  1]  I4)  C”1] 

a  push- forward  operation  yields  the  Cauchy  stress: 


a  =  — pl+ 


2  f 

1  f 


(Ci  +  C2/i)  B  —  C2B  +  C3  [exp  (/4  —  l)  —  l]  a  & 


a 


-1/3  (C1/1  +  2C2/2  +  C3  [exp  {U  -  1)  -  1]  h)  I] 


(51) 


5.3.2  Validation  of  Numerical  Method 

To  insure  that  the  stress  update  was  properly  implemented  into  the  computational  framework, 
numerical  tests  were  performed  and  compared  to  an  analytical  solution.  Single  element 
computations  of  uniaxial,  strip  biaxial,  and  equibiaxial  tests  were  performed  to  validate  the 
update  of  the  stress  for  the  implementation  against  theoretical  solutions.  If  the  stresses  are 
applied  in  the  x-y  plane  along  the  principal  axes  so  the  ozz  component  is  zero  and  the 
incompressibility  constraint  det  F  =  1  is  applied,  the  deformation  gradient  can  be  expressed  as 


F  - 


Ax  0  0 

0  0 


(52) 
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Then  using  equation  5 1 ,  the  stresses  are 


tfxx  =  2  [(Ci  +  C2h)  (Bxx  -  Bzz)  -  C2  (B2XX  -  B2zz)  + 
Cs  [exp  (74  -  1)  -  1]  74  {a2x  -  al)] 

(Tyy  =  2  [(<*  +  C2Il)  (Byy  ~  B „)  ~  C2  -  B^)  + 

C3  [exp  (74  -  1)  -  1]  74  (a2y  -  a*)] 


In  the  case  of  uniaxial  extension,  a}7  =  czz  =  0  and  'ky  =  'kz.  For  a  strip  biaxial  test,  azz  =  0  and  Xz  = 
1.  For  an  equibiaxial  test,  czz  =  0  and  Xx  =  Xy.  In  addition,  the  same  material  constants  were  used 
for  all  three  cases  including,  Ci  =  10.0,  C2  =  10.0,  C3  =  100.0,  and  k  =  10.0e9.  As  seen  in  figure 
18,  good  agreement  between  the  theoretical  and  finite  element  solutions  for  all  three  cases  of 
loading  is  achieved.  The  convergence  of  the  finite  element  solution  to  the  theoretical  solution 
demonstrates  that  the  stress  update  equations  were  implemented  correctly. 
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Figure  18.  Comparison  between  the  finite  element  results  and  analytic  solutions  for  three  different  loading 
conditions:  (a)  uniaxial,  (b)  strip  biaxial,  and  (c)  equibiaxial. 

5.3.3  An  Algorithm  for  Linking  DTI  and  Continuum  Mechanics 

The  transversely  isotropic  hyperelastic  model  previously  implemented  depends  on  the  fiber 
reinforcement  in  the  undeformed  configuration,  ao.  The  algorithm  developed  here  describes  how 
to  assign  each  finite  element  an  orientation  based  on  the  fiber  tractography  obtained  from  DTI,  as 
introduced  in  section  5.3.1.  This  is  complex  because  of  the  nonlinear  tractography  and 
unstructured  finite  element  shapes  within  the  three-dimensional  framework.  Algorithms  for  both 
hexahedral  and  tetrahedral  elements  are  described. 

First,  the  algorithm  divides  each  fiber  line,/i  (where  i  goes  between  0  and  the  number  of 
tractography  fibers),  into  individual  fiber  segments,  sj  (where  j  goes  between  1  and  the  number 
of  segments  within  a  given  fiber  line).  Each  fiber  segment  consists  of  two  fiber  segment  points 
with  coordinates  pi  and  P2:  s  =  P2  -  pi.  The  first  discretized  point  within  the  fiber  segment  is 
assigned  a  direction  based  on  the  vector  connecting  it  to  the  second  point.  If  the  fiber  point  exists 
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at  the  end  of  a  fiber  line,  that  terminating  point  receives  the  same  orientation  as  the  previous 
point.  A  schematic  of  the  discretization  process  is  shown  in  figure  19. 


Figure  19.  Fiber  discretization  and  notation  definitions  used  throughout  the  report. 

Next,  the  algorithm  determines  if  a  fiber  segment  within  a  given  tractography  fiber  overlaps 
spatially  with  a  finite  element.  In  order  to  determine  this,  various  cases  need  to  be  analyzed.  We 
begin  our  search  algorithm  at  the  fiber  segment  end  points.  For  each  element  in  the  volumetric 
mesh,  every  point  within  each  fiber  is  analyzed  for  the  cases  schematically  shown  in  figure  20. 
Case  1  exists  if  the  first  point  in  a  fiber  segment  is  contained  in  a  given  element.  Case  2  exists  if 
the  second  point  of  a  fiber  segment  is  contained  in  a  given  element.  Case  3  exists  if  both  end 
points  of  the  fiber  segment  are  contained  within  a  given  element.  Case  4  exists  if  a  fiber  segment 
intersects  the  faces  of  the  element,  but  its  fiber  points  are  not  contained  within  a  given  element. 
Case  5  exists  if  the  element  contains  none  of  the  fiber  points  and  the  segment  does  not  intersect 
the  element  faces.  Note,  that  there  is  an  unlikely  subset  of  case3  when  a  given  fiber  segment  lies 
parallel  with  a  single  element  face,  so  this  orientation  would  be  shared  by  any  adjacent  element 
also. 

For  all  of  the  cases  shown  in  figure  20,  determining  if  a  fiber  point  is  contained  within  a  finite 
element  is  complicated  by  the  unstructured  mesh  shapes.  For  example,  in  elements  with  adjacent 
faces  far  from  orthogonal  with  each  other,  a  fiber  point  may  lie  outside  of  the  faces  of  the 
element  even  if  the  fiber  point  spatial  coordinates  are  in  the  bounds  of  the  maximum  and 
minimum  element  nodal  coordinates.  Therefore,  an  algorithm  to  determine  if  the  fiber  point  is 
inside  of  the  element  is  implemented.  The  algorithm  takes  into  consideration  the  angles  of  the 
element  faces  explicitly  by  computing  normals  to  the  element  faces  and  the  coordinates  of  the 
fiber  point,  as  schematically  shown  in  figure  21. 
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Figure  20.  Schematic  of  the  various  cases  that  are  possible  when  determining 

if  a  fiber  segment  (within  a  given  tractography  fiber)  overlaps  with  a 
finite  element  in  space. 
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Figure  21 .  Method  used  to  test  if  a  specified  point  is  within  the  element.  Algorithms 

for  both  tetrahedron  and  hexahedron  elements  are  designed  and  implemented. 

A  slightly  more  complicated  algorithm  is  used  for  case  4  where  the  line  segment  intersects  the 
element  faces,  but  no  points  in  the  fiber  segment  reside  within  the  element. 
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As  schematically  shown  in  figure  22,  first  the  algorithm  calculates  the  normal  of  the  first  face,  n 
by  taking  the  cross  product  of  two  vectors,  which  share  the  same  vertex,  Pi. 


n  =  (P2  -  Pi)  x  (P3  -  P,) 


(55) 


Figure  22.  Method  used  to  test  if  a  line  intersects  the  faces  of  an  element  for  the  case  in  which  the  line  segment 
intersects  the  element  faces,  but  no  points  in  the  segment  reside  within  the  element  (see  Case  4  in 
figure  20). 

Next,  the  algorithm  takes  the  dot  product  of  the  distance  between  an  endpoint,  Li  or  L2,  and  the 
shared  vertex  Pi.  This  results  in  two  dot  products,  the  distances  from  Li  and  L2  from  Pi: 


X\  =  (Li  -  Pi)  •  (n) 

(56) 

z2  =  (L2  -  P0  •  (n) 

(57) 

If  xi  x2  >  0,  the  line  segment  never  intersects  the  plane  of  the  face.  If  xi  =  x2,  the  line  segment 
and  the  plane  are  parallel,  and  the  line  segment  can  be  disregarded  since  it  does  not  intersect  two 
different  faces  of  the  element.  If  the  two  preceding  conditions  do  not  exist,  the  code  will  find  the 
location  of  the  point,  P,  on  the  given  line  segment  which  intersects  the  plane. 


P  =  Li  +  (L2-Li)(-^~) 

X2  —  Xi 


(58) 
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The  point  P  is  the  intersection  point  that  exists  on  the  plane  of  the  element  face;  therefore,  the 
algorithm  must  test  if  P  also  lies  within  the  surface  area  of  an  element  face  (see  figure  22).  A 
point  is  inside  of  a  polygon  if  it  is  always  on  the  same  side  of  all  the  segments  that  create  a 
continuous  loop  around  the  edges.  For  a  given  fiber  segment  that  solely  intersects  the  element 
faces,  there  should  be  exactly  two  intersection  points  on  a  single  element. 

If  there  exist  elements  for  a  given  finite  element  mesh  where  no  fiber  tracks  overlap,  the 
elements  are  assigned  an  orientation  vector  of  zero,  ao  =  0.  Therefore,  the  transversely  isotropic 
material  condenses  to  an  isotropic  material. 

5.3.4  A  Scheme  for  Averaging  Multiple  Fibers 

The  previous  discussion  described  the  methods  used  to  find  the  existence  and  direction  of  a 
single  fiber  in  a  hexahedral  or  tetrahedral  element.  The  fiber  direction  assigned  to  an  element 
with  only  one  fiber  segment  is  given  by 

"°  “  INI  (51 

where  s  represents  the  orientation  vector  of  the  fiber  segment.  Depending  on  the  fiber  or  mesh 
densities,  there  could  be  multiple  fibers,  multiple  fiber  segments,  or  both  existing  in  a  single 
element,  as  schematically  depicted  in  figure  23.  The  current  effort  develops  an  average  scheme 
for  each  possibility  since  the  current  implementation  only  permits  one  orientation  vector  per 
element.  In  future  work,  we  hope  to  extend  the  model  to  have  the  capability  to  include  multiple 
fiber  orientations  within  a  single  finite  element,  which  will  require  additional  modifications  to 
the  theoretical  basis.  The  method  of  diffusion  spectrum  imaging  (DSI)  is  an  example  of  where 
this  functionality  would  be  useful  to  have,  since  DSI  is  capable  of  describing  fiber  crossings. 
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Figure  23.  Depending  on  the  fiber  or  mesh  densities,  there  could  be  multiple  fibers,  or  multiple  fiber 
segments,  or  both  existing  in  a  single  element  which  requires  an  averaging  scheme. 


In  the  case  of  a  single  fiber  with  multiple  fiber  segments  within  the  same  element,  an 
approximation  of  an  “average”  vector  direction  is  taken  because  the  material  model  uses  only 
one  fiber  direction  per  element.  By  simply  adding  the  vectors,  you  obtain  a  resultant  “average” 
direction,  which  provides  an  estimate  of  an  overall  fiber  orientation,  as  schematically  depicted  in 
figure  24. 


Figure  24.  The  sum  of  the  vectors  will  give 
an  average  orientation. 
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For  two  or  more  fibers  in  a  single  element,  the  assigned  orientation,  ao,  is  calculated  by  the 
following: 


max 

max  segments 
fibers  in  fiber  i 

ae  =  X  X  Sbll*l 


i=Q  j=l 

(60) 

a* 

a“=  IKII 

(61) 

where  Sj  are  the  orientation  vectors  of  the  fiber  segments  within  the  element  and  ae  represents  the 
average  of  Sj — two  averages  are  taken  with  each  element.  First,  the  vectors  from  the  segments,  Sj 
,  within  the  same  fiber  line  and  element  are  averaged.  Second,  the  average  of  all  fibers  within 
that  element  are  calculated.  Because  the  fiber  average  scheme  depends  on  the  element  size, 
future  efforts  will  explore  the  effects  of  finite  element  mesh  density.  A  course  finite  mesh 
decreases  computational  cost,  but  may  also  smear  out  the  complex  fiber  orientations.  However, 
there  are  regions  of  the  fiber  tractography  that  change  slowly  over  spatial  distance  so  a  course 
mesh  may  be  more  appropriate.  The  preceding  procedures  for  determining  an  element  orientation 
for  the  transverse  isotropic  model  can  be  summarized  by  the  global  algorithm  shown  in  figure 
25. 


for  cacli  element 
for  each  liber 

for  each  segment  in  fiber 

Determine  ease  for  fiber  points 
Store  segment  orientation  vector  s 
end 

Sum  segment  vectors  for  given  element  and  fiber 

end 

Compute  aL.  (Sum  orientation  vectors  for  all  fibers  in  element) 
Compute  a0 


Figure  25.  The  global  algorithm  used  to  obtain  final  vector  orientations  for  each  element. 

Although  this  method  will  automatically  take  into  account  the  vector  magnitudes,  one  limitation 
of  this  approach  involves  the  loss  of  orientation  information  due  to  the  averaging  scheme  for 
adjacent  fiber  segments  that  are  close  to  perpendicular.  For  example,  as  shown  in  figure  26,  if 
one  point  of  the  fiber  segment,  pi,  is  entirely  outside  of  the  element,  while  the  other  point  of  the 
segment,  p2,  is  inside  the  element  but  relatively  close  to  the  element’s  surface,  that  fiber  segment 
direction  will  receive  more  weight  than  it  should.  It  receives  the  same  amount  of  weight  as  it 
would  if  it  were  fully  contained  within  the  element.  Although  the  accuracy  of  the  fiber 
representation  is  limited  by  the  mesh  resolution,  it  is  possible  in  the  future  to  create  a  more 
accurate  averaging  of  fiber  directions  in  each  element. 
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Figure  26.  Averaging  fibers  that  are  not 

fully  contained  in  the  element  can 
create  a  misrepresented  average 
direction  for  the  individual  element. 


To  visualize  the  assigned  element  fiber  directions,  the  unit  vector  was  scaled  by  a  characteristic 
length  of  the  element  (diagonal  or  side  length)  and  added  to  the  center  coordinates  of  each 
element.  The  Visualization  ToolKit  (VTK)  output  file  describes  the  lines  connecting  the  center 
positions  and  endpoints  of  the  scaled  unit  vectors,  which  can  be  opened  up  in  Paraview  (Kitware) 
and  then  compared  against  the  original  dataset.  The  visualization  lines  are  not  meant  to 
completely  coincide  with  the  input  fiber  lines  but  rather  represent  the  directions  assigned  to  each 
element.  Since  the  visualization  lines  are  positioned  at  the  center  of  each  element,  the  lines  will 
be  in  close  proximity  to  the  fiber  lines  they  represent,  but  will  not  be  exactly  overlapping.  A 
higher  resolution  mesh  will  produce  a  more  accurate  representation  of  fiber  direction  assignment 
to  the  elements. 

5.3.5  Application  to  Modeling  White  Matter  Tissue  of  the  Human  Brain 

The  primary  objective  of  implementing  the  transverse  isotropic  model  was  to  incorporate  fiber 
directions  into  the  finite  element  models  for  application  to  anisotropic  biological  tissue.  Multiple 
scenarios  were  created  to  test  the  algorithms;  however,  instead  of  presenting  the  exhaustive  list 
here,  only  the  brain  model  is  shown.  Figure  27a  shows  the  input  fiber  tractography  superimposed 
with  a  finite  element  mesh  of  the  human  brain.  Figure  27b  shows  the  fiber  orientations  assigned 
to  the  elements,  obtained  from  the  fiber  tractography.  While  this  example  shows  qualitative 
agreement,  it  is  important  to  note  that  the  finite  element  mesh  and  the  fiber  tractography  are  not 
from  the  same  individual  at  this  time.  There  is  currently  a  major  focus  to  develop  models  from  a 
single  individual. 
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Figure  27.  DTI  input  data  versus  the  assigned  orientation  output. 

5.4  Conclusions 

The  work  presented  here  provides  a  novel  way  to  use  state-of-the-art  medical  imaging  to  inform 
a  three-dimensional  transversely  isotropic  model  of  biological  tissue.  After  testing  the  numerical 
model’s  implementation,  an  algorithm  was  created  in  order  to  read  fiber  data  and  output  the 
undeformed  fiber  directions  into  the  material  model.  The  initial  results  are  promising  since  they 
qualitatively  show  adequate  correlation  between  the  original  DTI  data  and  the  averaged 
directions  that  were  assigned  to  the  elements.  The  DWI-informed  modeling  can  be  extended  to 
other  biological  tissues  and  different  material  models.  Aside  from  using  DTI  for  imaging  the 
white  matter  of  the  brain,  it  can  also  measure  skeletal  muscle  fiber  directions  in  fixated  skeletal 
muscles.  Additionally,  the  model  could  be  useful  for  finite  element  modeling  of  non-biological 
fibrous  materials,  such  as  Kevlar  and  various  synthetic  polymers. 

As  mentioned  earlier,  future  work  will  include  finite  element  models  that  correspond  to  the  MRI 
data  and  DTI  data  derived  from  the  same  individual.  The  algorithm  could  be  improved  by 
constantly  updating  the  fiber  directions,  which  may  change  during  deformation,  although  this 
would  most  likely  increase  the  overall  time  required  to  run  the  simulation  and  may  not  be 
necessary  for  simulations  of  short  events,  such  as  blast  or  impact  loading.  Furthermore,  DSI 
images  maybe  incorporated  into  the  computational  framework  so  that  fiber  crossings  may  also  be 
included. 
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6.  Conclusion 


This  DSI  research  effort  in  Brain  Structure-Function  Couplings  aims  to  develop  a 
multidisciplinary,  multiscale  understanding  of  the  relationship  between  the  brain's  physical 
structure,  its  dynamic  electrochemical  functioning,  and  human  behavior.  This  report  provides  an 
overview  of  the  program  to  capture  the  general  direction  of  this  evolving  program  as  well  as  four 
project  descriptions  to  highlight  a  subset  of  Year  1  accomplishments  in  the  program. 

The  descriptions  for  Project  1  and  Project  2  highlight  the  development  and  assessment  of  a  time- 
evolving  functional  measure,  WPL1.  The  complementary  results  suggest  that  WPLf  may  have 
advantages  over  traditional  approaches  when  the  data  set  has  large-scale  movement  artifacts,  and 
time-evolving  WPLf  holds  promise  for  understanding  brain  function  in  complex  operational 
environments  with  real-world  eye  and  body  movements,  providing  a  critical  tool  for  the 
development  of  neurotechnologies  that  could  monitor  neural  data  in  real  time  to  augment 
Soldier-system  performance. 

The  descriptions  for  Project  3  and  4  highlight  accomplishments  under  the  program’s  two 
modeling  efforts.  The  electrochemical  structure-function  modeling  results  indicate  that  the  three 
phase  measures  tested  (ImC,  dWPLl,  and  PLV)  do  capture  the  temporal  dynamics  of  the 
communication  between  nodes,  but  even  in  two-node  networks,  the  results  indicate  that  network 
structure  strongly  influences  the  measure  results  and  suggests  relationships  between  structure 
and  function  need  to  be  explored  further.  The  biomechanical  structural  modeling  effort  results 
describe  a  method  for  incorporating  empirically  collected  white  matter  tract  data  (DW1)  into  a 
finite  element  model.  This  work  lays  the  critical  foundation  needed  to  increase  the  fidelity  of  the 
brain  tissue  response  in  simulations  of  tissue  damage  due  to  blast  or  blunt  forces. 

Collectively,  these  four  projects  capture  the  research  enhancements  derived  through  cross- 
Directorate  collaboration.  The  simulated  EEG  data  can  be  used  to  validate  and  better  understand 
functional  measures  that  are  tested  on  empirical  datasets.  The  imaged  brain  structure  is 
incorporated  into  a  finite  element  model  to  increase  the  fidelity  of  the  simulated  tissue  response 
properties.  Critically,  in  our  first  year,  each  of  the  four  efforts  has  successfully  incorporated  a 
graph  theoretic  framework  to  ensure  advancements  between  areas  can  be  integrated  as  the 
research  develops  in  the  years  to  come. 
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ARL 

BSFC 

CISD 

CTA 

DRC 

DSI 

DTI 

DWI 

dWPLI 

EEG 

FA 

FIR 

fMRI 

fNIR 

HRED 

Hz 

ICA 

ICB 

ICs 

ImC 

LFP 

MEG 

MRI 

MRMC 


U.S.  Army  Research  Laboratory 

brain  structure-function  couplings 

Computational  and  Information  Sciences  Directorate 

Collaborative  Technology  Alliances 

Dynamics  Research  Corporation 

Director's  Strategic  Initiative 

diffusion  tensor  magnetic  resonance  medical  imaging 

diffusion-weighted  imaging 

debiased  weighted  phase  lag  index 

electroencephalography 

fractional  anisotropy 

finite  impulse  response 

functional  magnetic  resonance  imaging 

functional  near-infrared  spectroscopy 

Human  Research  and  Engineering  Directorate 

hertz 

independent  component  analysis 
Institute  for  Collaborative  Biotechnologies 
independent  components 
imaginary  component  of  coherence 
local  field  potential 
magnetoencephalography 
magnetic  resonance  imaging 
Medical  Research  and  Materiel  Command 
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MURI 


Multidisciplinary  University  Research  Initiative 
NFFT  nonequispaced  fast  Fourier  transform 

NMM  neural  mass  model 

OSTP  Office  of  Science  and  Technology  Policy 

PLI  phase  lag  index 

PLV  phase  lag  value 

PSD  power  spectral  density 

UCSB  University  of  California  at  Santa  Barbara 

USMA  United  States  Military  Academy 
VTK  visualization  toolkit 

WMRD  Weapons  and  Materials  Research  Directorate 
WPLI  weighted  phase  lag  index 
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